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Introduction 

Simulation  of  analog  hard  faults  in  piecewise  linear  networks  could  be  done 
efficiently  using  the  N-port  theory  for  formulating  the  network  equilibrium  equations 
then  using  the  complementary  pivot  theory  to  solve  for  the  node  voltages  [Ij.  For  a 
network  of  n  ports,  a  set  of  n  hybrid  equations  could  be  obtained  using  Lin’s  method 
[2],  These  equations  act  as  a  data  base  for  simulating  different  faults  without  the  need 
to  reformulate  new  equilibrium  equations.  Every  time  a  fault  is  simulated,  only  a  sub¬ 
set  of  these  equations  have  to  be  solved  using  the  Lemke’s  complementary  pivot  algo¬ 
rithm  [3]  to  obtain  the  node  voltages  of  the  preselected  test  nodes.  An  analog  fault  dic¬ 
tionary  could  be  obtained  by  logical  means  whereby  faults  are  identified  by  numerical 
codes  [4].  A  FORTRAN  program  is  included  in  this  document  to  achieve  the  above 
tasks.  The  program  user  is  expected  to  provide  a  piecewise  linear  model  of  the  network 
to  be  analyzed. 

Hardfault  Modelling: 

Open  circuited  inductances  and  short  circuited  capacitances  are  conveniently  simu¬ 
lated  without  extra  work,  whereas  simulating  faults  in  other  elements  require  the  addi¬ 
tion  of  switches  which  are  normally  open  (NO)  or  normally  closed  (NC).  Nodes  which 
are  chosen  as  test  nodes  must  be  considered  to  form  a  zero  value  current  source 
directed  from  the  test  node  to  the  ground  node.  This  is  for  compatibility  with  the  for¬ 
mulation  to  make  it  possible  to  solve  for  the  node  voltages  by  letting  the  zero  value 
current  source  be  a  current  port  in  the  n-port  formulation.  The  basic  idea  is  to  pull  the 
following  elements  out  of  the  network  to  be  considered  as  current  or  voltage  ports  as 
3hown  in  Fig.  1. 

1.  independent  voltage  sources. 

2.  independent  current  sources. 

3.  normally  open  switches. 

4.  normally  closed  switches. 

5.  ideal  diodes. 

6.  zero  value  current  sources. 

7.  inductances. 

8.  capacitances. 

Fig.  2  shows  possible  ways  of  modelling  some  common  faults.  It  can  be  seen  that  inser¬ 
tion  of  a  normally  closed  switch  requires  the  addition  of  a  new  node.  Similar  ideas  can 
be  used  in  simulating  other  types  of  faults.  For  example,  figure  3  suggests  a  way  for 
simulating  an  open  circuited  base  or  a  short  circuited  base  emitter  juuction  in  a  piece- 
wise  linear  model  of  a  bipolar  transistor. 


Fig.  3.  Simulation  of  open  base  and  short  circuited  base  emitter 
junction  in  piecewise  linear  models  of  bipolar  transistor. 
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Description  and  Execution  of  the  program: 

The  program  requires  the  network  description  to  be  stored  in  a  file  named  data  in 
a  special  format.  A  program  “cctl”  has  been  written  to  accept  the  network  description 
from  the  user  interactively  and  prepare  the  required  file.  Any  changes  or  corrections  in 
the  data  file  afterwards  must  preserve  the  required  format.  The  program  “cctl”  is 
compatible  with  the  UNIX  f77  compiler.  To  prepare  “cctl”  from  the  source  file  cctl.f, 
use  the  following  command  sequence 

(77  cctl.f 

mv  a.out  cctl 

cctl 

The  current  version  of  the  main  program  is  compatible  with  the  MNF  compiler 
that  is  available  on  the  PVCC  CDC  6600  computer  system  under  the  control  of  a  dual 
MACE  operating  system.  Files  can  be  moved  back  and  forth  between  the  ECN  com¬ 
puters  and  the  PUCC  computer  using  the  UNIX  “pjs”  command.  For  complete  infor¬ 
mation  about  pjs,  see  the  help  file  available  on  the  UNIX  by  typing  Shelp  pjs.  Gen¬ 
erally  the  following  forms  of  the  command  are  enough  to  put  a  file  in  or  get  a  file  from 
the  PFILES  storage  in  PUCC: 

$pjs  -  put  -  use  [id]  filenamel:  filename2 
or 

$pjs  -  get  -  use  [id]  filenamel:  filename2 
where 

id:  is  the  user  id  in  the  PUCC. 
filenamel:  is  the  UNIX  file  name 
filename2:  is  the  PFILES  file  name 

If  the  filename  where  data  is  to  be  transferred  is  not  in  the  system,  a  new  file  will 
be  created  and  given  the  specified  name.  The  computer  will  then  prompt  by  asking 
about  the  account  number  and  password  on  the  PUCC  computer.  The  main  program 
IIAFDIC  uses  the  following  4  main  subroutines  that  are  available  in  4  separate  files: 

1.  HYBRID 

2.  LEMKE 

3.  AMBSET 

4.  FTCODE 

None  of  these  subroutines  requires  simultaneous  availability  of  another  subroutine  in 
the  core  which  makes  possible  overlay  loading  in  a  smaller  core  size  if  so  desired. 
Without  overlay  loading,  the  required  core  size  is  134000  words. 


The  simplest  way  of  execution  given  the  source  files  is  to  compile  every  file 
separately  and  keep  a  copy  of  the  binary  object  files  resulting  from  compilation  to  be 
loaded  individually  any  time  a  program  is  needed. 

The  deck  required  for  compiling  HYBRID  is: 

12345,  ABC,  MF 77000,  CM77000. 

PFILES  „  HYBRID. 

RFL  (77000) 

MNF  (I  =  HYBRID,  B  =  BHYBRD,  N,  L  =  0 
PFILES,  PUT,  BHYBRD. 

#EOR 

#EOF 

The  parameter  N  in  the  MNF  command  will  prevent  execution  and  L  =  0  will  suppress 
listing.  The  maximum  field  length  and  central  memory  size  required  for  compiling  any 
subroutine  are  those  required  for  HYBRID.  Therefore  the  same  deck  can  be  used  for 
compiling  other  subroutines.  The  names  of  the  binary  files  available  now  on  the  system 
are  BHAFDIC,  BHYBRD,  BLEMKE,  BAMBST,  BFTCOD.  It  has  to  be  noted  that  the 
commands  PFILES  and  MNF  automatically  set  the  field  length  to  the  default  values  of 
15000  and  45000  respectively.  The  45000  words  field  length  is  not  enough  for  compiling 
some  subroutines.  Therefore  it  is  necessary  to  use  the  RFL  command  for  adjusting  the 
field  length.  The  deck  required  for  that  is: 

ABC,  12345,  MF 60000,  CM60000. 

PFILES,,  AMBSET. 

RFL, 60000. 

MNF  (I  =  AMBSET,  B  =  BAMBST,  N,L  =  0) 

PFILES,  PUT,  BAMBST. 

PFILES  „  FTCODE. 

RFL,  60000. 

MNF(I  =  FTCODE,  B  =  BFTCOD,  N,  L  =  0) 

PFILES,  PUT,  BFTCOD. 

PFILES,,  LEMKE. 

RFL,  60000. 

MNF  (I  =  LEMKE,  B  =  BLEMKE,  N,L  =  0) 

PFILES,  PJT,  BLEMKE. 

PFILES,,  TESTN. 

RFL,  60000. 

MNF  (I  =  TESTN,  B  =  BTESTN,  N,L  =  0) 

PFILES,  PUT,  BTESTN. 
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The  deck  required  for  execution  is 

12345,  ABC,  MF 134000,  CM134000. 

PFILES  „  BTESTN. 

PFILES  „  BHYBRD. 

PFILES  „  BLEMKE. 

PFILES  „  BAMBST. 

PFILES  „  BFTCOD. 

RFL,  134000. 

LOAD,  BHYBRD. 

LOAD,  BLEMKE, 

LOAD,  BAMBST. 

LOAD,  BFTCOD. 

LOADX,  BTESTN. 

REWIND,  OUTPUT 

PFILES  (PUT,  RESULT,  X  =  OUTPUT) 

#EOR 

#EOF 

The  last  two  commands  will  store  the  output  in  a  file  called  RESULT  which  can  be 
printed  on  the  line  printer  afterwards  at  the  user’s  convenience,  or  kept  for  inspection 
in  the  PFILES  storage. 

Example 

The  video  amplifier  in  Fig.  4  has  the  piecewise  linear  model  shown  in  Fig.  5,  where 
branch  numbers  are  enclosed  in  circles  while  node  numbers  are  written  beside  the 
corresponding  nodes.  The  user-computer  dialogue  as  controlled  by  “cctl”  is  shown 
next.  Next  to  it  is  the  data  file  produced.  Note  the  second  line  in  the  file  which  con¬ 
tains  14  integer  numbers  and  a  floating  point  number.  This  line  gives  the  numbers  of 
branches  in  the  14  permissible  types  of  two  terminal  elements  followed  by  the  ambi¬ 
guity  voltage  range.  If  the  user  decides  to  add  or  delete  any  branches  before  running 
the  program,  he  has  to  modify  the  number  of  branches  accordingly.  The  numbers  in 
the  second  line  follow  the  same  sequence  which  appears  in  the  question  “branch 
type?(e,L,..etc.)”.  The  program  output  follows  the  listing  of  the  data  file. 

Simulating  Multiple  Faults 

If  several  branches  are  to  be  considered  simultaneously  faulty,  they  should  be  all 
assigned  the  same  fault  number  as  fault  0  in  the  example.  If  the  fault  in  some  branch 
is  not  required  to  be  simulated  at  all,  the  field  corresponding  to  the  fault  number 
should  be  skipped  when  replying  to  the  question  concerning  this  branch,  as  in  branch 
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48  in  the  example. 

Branch  Polarities: 

The  end  nodes  of  all  branches,  except  diodes  and  controlled  branches,  can  be 
entered  in  either  order.  The  values  of  voltage  and  current  will  be  adjusted  accordingly 
to  be  either  positive  or  negative. 

For  diodes,  the  node  connected  to  the  n  terminal  should  be  entered  first  to  make 
both  current  and  voltage  positive  in  either  case  of  being  on  or  off,  otherwise  he  comple¬ 
mentary  problem  will  not  be  feasible.  The  program  assumes  the  port  voltage  and 
current  convention  as  shown  beside  the  diode  element  in  figure  1. 

For  control  branches,  the  order  of  entering  the  end  nodes  must  be  consistent  with 
the  order  of  entering  the  end  nodes  of  the  controlled  source  according  to  the  current 
and  voltage  conventional  directions  shown  in  figure  1. 

Numbering  Nodes  and  Branches: 

Branches  can  be  given  any  numbers  which  are  not  necessarily  consecutive.  How¬ 
ever  node  numbers  must  be  consecutive  starting  from  zero,  which  is  to  be  the  ground 
node.  Branch  names  can  be  up  to  4  characters  which  must  start  by  one  of  the  14  acro¬ 
nyms  (e,j,...)  shown  in  the  example.  The  names  of  the  diodes  and  the  test  node  parts 
are  the  ones  to  be  used  in  the  program  output. 

References: 
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problem  title! 

a  new  video  amplifier  circuit 
branch  type(e*J*r*g,cc*cv*vc*vv*vt*L*c*nc*no*d>: 

e 

branch  no* f rom, to, bat tery  value 

7,6»0*-28 

branch  type<e*]*r*g*cc*cv*vc*vv*vt*L*c*nc*no*d): 

e 

branch  no  * f rom « to « ba t ter y  value 
14*28*0 .28 

branch  t  ype(e*  j  *r«g«cc  *cu.vc«vv»vt  *l*c*nc«no«(j)  : 
e 

branch  no* f rom* t o* bat t ery  value 
5  4  5  .7 

branch  tyDe(e«j«rfg*cc*cv«vc*vv*vt*L*c*nc*no*d): 
19*11, 12*. 7 

branch  type<e*j*r*q«cc.cv*vc.vv*vt*L«c*nc*no*d>: 
e 

branch  no* f rom* t o* ba t t er y  value 
34. 17, 18, .7 

branch  type(e*j*r*g*cc«cv*vc*vv*vt*L*c*nc*no.d): 
e 

branch  no* f rom* t o» ha t t ery  value 
51*24, 25, .7 

branch  tyoe(e*j*r*g*cc*cv*vc*vv*vt*L«c*nc*no*d>: 
ol 

branch  no.frorr.to 
4,4,3 

branch  type(e«]*r*g*cc*cv*vc*vv*vt ,L*c*nc*no*d): 
d  2 

branch  no.from.to 
10*7,5 

branch  type(e*  j,r,g,cc*cv*vc*vv*vt*L*c*nc*no*o): 
d3 

branc  h  no*  f  rom, t  o 
18*11,10 

branch  tyoe(e*]*r*g*cc*cv*vc*vv*vt*L*c*nc*no*d>: 
d4 

branch  no.from.to 
25,14,12 

branch  tyoe(e*]*r*g*cc*cv*vc*vv*vt*L*c*nc*no*d)  : 
d5 

branch  no.from.to 
33*17*16 

branch  type(e*j*r*g*cc«cv*vc*vv*vt*L*c*nc*no*d>: 
d6 

branch  no»from»to 
42*20*18 

branch  type(e*]*r«g*cc*cv*vc*vv*vt*L*c*nc*no*d): 
d7 

branch  no.from.to 
50,24,23 

branch  typete*  j*r*g*cc*cv*vc*w*vt*L*c*nc*no*d): 
d8 

branch  no.from.to 
58*27,25 

branch  tyoe(e*j*r*c*cc*cv*vc*vv*vt*L*c*nc*no*d>: 
nc 

branch  no* f r om, t o* f au l t  no(1f  o.c  Is  to  be  simulated) 
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1. 1.2,2 

branch  t  ype( e « ],r,gtcc»cvtvc.vv,vt,L,c,nc,no,o): 
nc 

branch  no « f rom « t o * f aul t  no(1f  o.c  Is  to  oe  simulated) 

16*7,9,3 

branch  type<e,J,r»g»cc»cv»vc»vv*vt»L,c»nc*no»d>: 
nc 

branch  no, f r om * t o» f au  1 1  no(1f  o.c  Is  to  be  simulated) 

26. 14.8.6 

branch  type(e,],r,g,cc,cv,vc,vv,vt,L,c,nc,no,d): 
nc 

branch  no ,f rom , t o , f aul t  no(1f  o.c  Is  to  be  simulated) 
31,15,29,4 

branch  type(e, j*r,g,cc«cv«vc,vv,vt*L»c,nc,no»d): 
nc 

branch  no , f r om , t o , f aul t  no(1f  o.c  Is  to  be  simulated) 

48,22,30, , 

branch  type(e,]«r,g*cc,cv,vc,vv,vt,L,c,nc,no,d): 
nc 

branch  no, from, to, fault  no(1f  o.c  Is  to  be  simulated) 

60.27.8.7 

branch  type(e,j,r,g,cc«cv«vc,vv*vt*L«c,nc«no,d): 
r 

branch  no, f rom , t o, res  1 s t ance  value 
61,1,0,1200 

branch  type(e,],r,g,cc,cv,vc,vv,vt,L,c,nc,no,d): 
r 

branch  no , f rom , t o  ,  res  1  stance  value 
3,2,3,500 

branch  typete, j,r»g*cc*cv*vc,vv»vt,L,c,nc»no»d)! 
r 

branch  no,  f  roir ,  t  o,  re s  1  s  t  ance  value 
6,5,6,5670 

branch  tyoe(e,],r,g,cc,cv,vc,vv,vt,L,c,nc,no,d>: 
r 

branch  no, f rom, t o , res  1 s t an ce  value 
12,7,8,3000 

branch  type(e, j,r,g,cc,cv,vc,vv,vt,l,c,nc,no,d>: 
r 

branch  no, f rom, to, res Istance  value 
17,9,10,1000 

branch  tyoe(e,J,r,g,cc,cv,vc,vv,vt,L,c,nc,no,d): 
r 

branch  no, f rom, to , resl s t ance  value 

20,12,0,1200 

branch  type(e«  j,r,g,cc,cv,vc,vv,vt,L,c,nc,no,d): 
r 

branch  no, from, to, res Istance  value 
21,13,0,50 

branch  type(e,],r,g,cc,cv,vc,vv,vt,L,c,nc,no,d): 
r 

branch  no, from, to, r^s  Istance  value 
13,28,a, 78 

branch  type(e,],r,g,cc,cv,vc,vv,vt,L,c,nc,no,d): 
r 

branch  no, from, to, resistance  value 
28,12,15,330 

branch  type(e,],r,g,cc,cv,vc,vv,vt,L,c,nc,no,d>: 
r 

branch  no, f r om , t o, r f s 1 s t ance  value 
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39,12,22,330 

branch  type(e,J,r»g,cc»cv»vc,vv»vt»l,c»nc»no,d)! 

r 

branch  no, f row, to, res fstance  value 
32*29,16,1000 

branch  type(e»j»r»g»cc»cv»vc»vv*vt*L»c»nc»no»d>: 
r 

branch  no , f rom , t o , r es 1 s t ance  value 
35*18,19,500 

branch  type(e«]*r,g,cc,cv,vc*vv.vt«L«c,nc,no,d): 
r 

branch  no , f rom , t c , re s 1 s t anc e  value 
37,19,0,1 200 

branch  tyoe(e,J,r,g,cc,cv,vc,vv,vt,L,c,nc,no,d): 
r 

branch  no , from , t o , res  1 st ance  value 

90.20.8.1000 

branch  type(e,J«r,q,cc,cv,vc»vv»vt»L»c,nc»no»d): 
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2.264e*01 


000000000001000 


7 

15 

2.425e*01 

2 .490e*01 

8 

16 

2.084e*01 

2. 164e*01 

node 

no: 

v  1 1 2 

set 

center  from 

t  0 

1 

1 

1.023e*01 

1. 083e*01 

2 

2 

2.271e«-01 

2.351e«-01 

3 

3 

-4.000e-01 

4 •  000  e-01 

4 

6 

2 • 39 1 e*  0  0 

3.000e>00 

5 

4 

1.083e*01 

1.143e*01 

6 

10 

3.415e«>00 

4.022e*00 

7 

13 

3.000e*00 

3.415e*00 

8 

14 

9.425e*00 

1. 023e<-01 

9 

15 

8.420e*C0 

9.220e*00 

10 

16 

5.395e*00 

6.195e*00 

node 

no : 

vt  1 3 

set 

center  from 

t  o 

1 

1 

-4.000e-01 

4 • 00  Oe-0 1 

2 

13 

2«808e*00 

3«608e+Q0 

node 

no : 

vt  15 

set 

center  from 

to 

1 

1 

1.020e«-01 

1.081e«-01 

2 

2 

2«160e+01 

2.240e*01 

3 

3 

-4»000e-01 

4 . 000e-01 

4 

6 

2.385e*00 

2.993e*00 

5 

4 

1.081e*01 

1  •  1 43e*0 1 

6 

10 

3.408e*00 

4«014e^00 

7 

11 

8.317e*G0 

9.1l7e»00 

000000000030010  -2 

000000000000001 

set  code 

100001010110000 

010000100000000 

001000000000000 

000100000000000 

OOOOIOOCOOOOOQQ 

000000001000000 

000000000001000 

000000000000100 

000000000000010 

000000000000001 

set  code 

111111111110111 

000000000001000 

set  code 

100001010010000 

010000100000000 

001000000000000 

000100000000003 

000010000000000 

000000001000000 

000000000100000 


8 

13 

2.993e«-00 

3.408e*00 

9 

14 

9.297e*00 

1. 010e  +  01 

10 

15 

6.747e«-00 

7.547e*00 

11 

16 

5.380e«-00 

6. 1 80e*00 

node 

no : 

vt  18 

set 

center  from 

t  0 

1 

1 

9.41 3e*00 

1  .  021e*01 

2 

2 

1.755e*01 

i.835e^01 

3 

3 

-4.00 Oe -01 

4. 000e-01 

4 

6 

1  •  722e>  00 

2. 27?e«.00 

5 

10 

2.684e*00 

3.289e*00 

6 

11 

8.31 7e*00 

9.117e*00 

7 

13 

2.273e*00 

2. 6R4e*00 

8 

15 

9.767e-01 

1 • 722e*00 

9 

16 

4.637e*00 

5. 437e*00 

node 

no : 

vt  19 

set 

center  from 

to 

1 

1 

6.540e*00 

7 • 32  7e*0  0 

2 

1  •  227e*0 1 

1 • 307e+Cl 

3 

3 

-4.000e-01 

4. OOOe-Ol 

4 

6 

1.21 5e >00 

1 .749e*00 

5 

10 

1  •  74  9e  ♦  0  0 

2.439e*00 

6 

11 

5.753e*00 

6.540e*00 

7 

15 

5.7 1 8e-0 1 

1.215e*00 

8 

16 

3. 155e*00 

3.955e*00 

000000000001030 
000000000000100 
00000000  "0-00.  0010 
000000000000001 


set  code 

100001010010000 

010000100000000 

001010000000000 

000100000000000 

000000001000000 

000000000100100 

000000000001000 

000000000000010 

000000000000001 


set  code 

100001010010000 

010000100000000 

001010000000100 

000100000001000 

000000001001300 

000000000100000 

000000000000010 

000000000000001 


node  not  vt20 


center 


from 


set  cod** 


1 

1 

2.016e*01 

2. 096e*01 

2 

2 

1.755e*01 

1 .835e*01 

3 

3 

2.723e*01 

2.8  0  3e*0 1 

4 

6 

2.581e*01 

2.63He*01 

5 

4 

2.638e*01 

2.707e*0l 

6 

10 

2.511 e* 01 

2.5«le*01 

7 

14 

8.20  7e  ♦  0  0 

9.007e*00 

8 

15 

9. 767 e- 01 

1.777e*00 

node 

no : 

vt  22 

set 

center  from 

t  0 

1 

1 

1  •  024e*  0 1 

1.  101e*01 

2 

2 

2.26Se*0 1 

2. 348e ♦ 0 1 

3 

3 

-4. Q00e-0 l 

4,0  0  0e “01 

4 

6 

2 . 388e  >00 

2. 996e*00 

5 

7 

9.47 3e *00 

1.024e*01 

6 

10 

3.411e*00 

4 • 018e*00 

7 

13 

2.996e*00 

3.411e*00 

8 

15 

8.408e*00 

9.208e*00 

9 

16 

5.035e*00 

5.835e*00 

node 

no : 

»/t25 

set 

center  from 

to 

1 

1 

9.469e*00 

1 • 027e*01 

2 

2 

2. 188e*01 

2 . 268e*0  1 

3 

3 

-4.000e-01 

4. 000e-01 

4 

6 

1.679e*00 

2 • 286  e*00 

5 

7 

6.640e*00 

7.440e*00 

6 

10 

2.699e*00 

3. 276e*00 

7 

13 

2 .286e*00 

2.699e*00 

100001000011000 

010000100000001 

001000000000000 

000100010100000 

000010000000000 

00000000100000C 

000000000000100 

000000000000010 

set  code 

100010010100000 
010000100000000 
QOIOOOOOOOOCOOO 
0001000000000' 
000001000010100 
00000  OCi'lOtfCOOO 
000000000001000 
00000000000001G 
000000000000001 

set  code 

100010010110000 

010000100000000 

001000000000000 

000100000000000 

000001G0000003U 

000000001000000 

000000000001000 


8 

14 

8.672e+00 

9.46=e*0C 

9 

1  5 

7.672e*00 

8«472e*01 

10 

16 

3.276e*00 

4.047e*00 

node 

no  1 

vt2o 

ooococcoooaoioo 

000000003000010 

000000000000001 


set  center  from  to 

1  1  -4.00 Oe -01  4.000e-01 

2  16  3.247e*00  4.047e«-00 


set  code 

111111111111110 

000000000000001 


fault 

fault 

code 

f  1 

vt5 

v  1 7 

v  1 8 

vtl2 

vt  1 3 

f  9 

1 

1 

1 

1 

1 

f  12 

1 

1 

1 

4 

1 

f  7 

1 

1 

5 

1 

1 

f  2 

1 

e 

5 

1 

1 

f  9 

2 

2 

2 

2 

2 

f  3 

2 

2 

2 

2 

3 

f  6 

3 

3 

3 

3 

1 

f  4 

4 

4 

4 

4 

1 

fio 

5 

1 

1 

5 

1 

f  11 

6 

6 

6 

6 

1 

f  13 

7 

1 

1 

4 

1 

f  1 4 

8 

7 

7 

1 

1 

f  15 

9 

8 

5 

7 

1 

f  16 

10 

9 

8 

8 

1 

#eor 

#eof 

11 

10 

9 

2 

1 

Appendix  I 

Copyright,  Purdue  Research  Foundation,  1083. 


program  cctl  for  accepting  network  topology  for  fault  simul- 
tion  and  producing  a  special  format  output  file  containing 
the  network  information.  The  output  file  name  is  specified 
by  the  user.  To  oe  compatible  with  the  "hafdic"  oroqram, 
tne  output  file  should  be  named  "data". 

internal  variables 

ne  number  of  Independent  voltage  sources 

nl  number  of  Inductances 

nnc  number  of  normally  closed  switches 

no  number  if  ideal  diodes 

nr  number  of  resistances 

nq  number  of  conductances 

ncc  number  of  current  controlled  current  sources 

ncv  number  of  current  controlled  voltaqe  sources 

nvc  number  of  voltage  controlled  current  sources 

nvv  number  of  voltaqe  controlled  voltaqe  sources 

nvt  number  of  test  nodes 

nno  number  of  normaly  open  switches 

nc  number  of  capacitances 

nj  number  of  independent  current  sources 

or  vector  containing  branch  numbers 

nft  vector  containing  fault  numbers 

bat  vector  containing  batterry  values 

cur  vector  containing  values  of  Independent  current  sources 

nfrom  vector  containing  source  nodes  of  the  cor  re soona 1 ng 
branches 

nto  vector  containing  destination  nodes  of  the  correspond¬ 

ing  branches 

icont  array  containing  control  branches. 

x  array  containing  values  of  resistancest conductances 

and  control  branches, 
title  array  containing  problem  title 

ch  array  containing  branh  names 

filnam  character  variable  containing  file  name  soecifleo  by 
the  user  for  storing  data 

limitations 

*0  nodes 
90  branches 

40  voltage  Independent  sources 
40  current  Independent  sources 
40  dependent  sources  (all  kinds) 

60  resistances  or  cinductances 


Integer  nf t < 90  )  tbr < 90 ) 

dimension  bat(30)«cur(30)«nfrom(90)»ntoC50)»1cont<4»90) 

dimension  x(f.*40) 

character *80  title 

character*4  1stoo*ch<90) 

cnaracter*20  filnam 

character*!  1e«15?11»1rt1qtil*1d*ic 

character*2  1nc»1vt»1no«1ect1cv*1vc*1vv 


data  icciicvtivciivv/'cc'i'cv't'vc'i'vv’/ 

data  1stop»te*1j»i1*1r*1q/'stop'»'e**'J'»'1***r»*'g'/ 

data  1l*1nc«1a«tvt»1no»1c/'L'*'nc'»'d'*'vt**'no'»*c'/ 

ne  =  0 
n  l  =  0 
nnc  =  0 
nd  =  0 
nr  =  0 
ng  =  0 
n  c  c  =  0 
new  =  0 
n  vc  =  0 
nvv  =  0 
nvt  =  0 
nno  =  0 
nc  =  0 
n )  =  0 
nD  =  0 

print*  *  'problem  title:* 
read( 5* • ( a80 ) * )  title 

do  5  1=1*90 
nf t (1  1  =  0 

do  100  k=l,9C 

pr1nt**'branch  tyoe(etJ«r«g*cc*cv*ve*vv*vttL»c»nc*notd) 
r ead( 5 • * ( a4 ) • )  ch(k> 

1 f ( ch ( k ) «eq. 1 s too )  qo  to  110 
If(chUXr.l).eq.ie)  then 
ne=ne*l 
np=no*l 

pr  Int *» 'branch  no « * r on* t o* ba 1 1 e r y  value' 
read**  br(k)»nfrom(k>*nto<k>»bat(ne> 
elself (ch(k)(i:i).eq.1l)  then 
nl=nl*l 
np=np*l 

print •* »b ranch  no* fromtto » fault  no(1f  o.c  Is  to  be  simulated) 

read**br<k)*nfrom(k)*nto(k)»nft<np> 

elsel f ( ch< k)< 1 :2> .eq.lnc)  then 

nnc  =nn  c ♦ 1 

no=np+ 1 

pr  1  nt *♦' branch  no  * f r om * t o ♦ f au  1 1  no(1f  o. c  Is  to  be  simulated) 

read**br<k)»nfrom(k)*nto<k>»nft(no> 

e l se 1 f < ch < k ) ( 1 : 1 ) . eq .  1 d )  then 

nd=nd* 1 

no=np*l 

pr 1 nt * * 'branc h  no«from*to* 

reaa**br(k)*nfrom(k)»nto(k) 

e Isel f ( ch ( k ) ( 1 : 1 ) • eq.  1  r )  then 

pr 1 nt *»'  branch  no « f rom » to  *  res  1 s t ance  value' 

nr  =  nr*  1 

reac**br(k)»nfrom(k)*nto(k)*x(l*nr) 

elsel  f (ch (k  )  C 1 :  1 ) .eq.  1q)  then 

print* » » oranch  no, from* to* conduct ance  value' 

ng  =ng  + 1 

reaa**br<k)*nfrom<k>*nto(k)»x(2*ng> 
elsel f < ch ( k ) ( 1 : 2 ) • eo.  Icc )  then 

pr Int *» 'branch  no* f rom»to*control  branCh*cc  value' 


nc  c  =  nc  c* 1 

read**br(k)*nfrom(k)*nto(k)»icont(l*ncc)»x(3«ncc> 
e l se i f ( c h ( k ) ( 1 : 2 > • eq •  i  c v >  then 

pr 1 nt * » * br anc h  no  * f rom * t o ♦ c on t r o l  branch*cv  value* 
ncv=ncv*l 

reaa*»  br(  k)  »  n  f  r  c  m  (  k  )  *  n  t  o  (  k)  «  Icont  (  2*ncv  )  *  x  (  4  •  n  c  v  ) 
e l se 1 f ( c h ( k ) ( 1 : 2 ) . e a .  i  vc  )  then 

or int* *' branch  no « f r om , t o * c on t r o l  brunch*vc  value* 
nv  c=nv  c*l 

read*tDr( k) *  n f  rom ( k ) *  nt  o ( k  ) ,  icont ( 3*nvc)*x(5*nvc) 
e l s e i f ( ch ( k ) (  1 :  2 > . e a  .  1 vv  )  then 

p r i n t * «  •  o r a nc h  no  *  1 r on 1 1 o * c on t r o  l  oranch*vv  value' 
nvv=nvv* 1 

read*  *  b  r  (  k ) *  nf  r  om ( k) *nt  o  ( k ) » i cont (4*nvv)*x(6*nvv> 
e 1 se i f ( ch ( k ) ( 1 : 2 ) • eq.  i v t >  then 
nvt  =nvt  *1 
np=np*l 

pr int *♦' branch  no»from*to' 
read*«br( k) *  n  f  rom( k ) »nto( k) 
elsei f ( ch ( k ) ( 1 :  2 ) .eo  •  i no  )  then 
nno  =  nno  *1 
np  =  np*  1 

p r i n t * « • b ranc h  no* from* to* fault  nolif  s.c  is  to  be  simulated)' 
read*»br(k)*nfrom(k) *nto(k) *nft(np) 

elseif(ch(k)(i:i).eo.ic*and.ch(k)(i:2).ne.icc.and.ch(k)(i:2).ne 
1  cv)  then 
nc  =  nc ♦ 1 
np  =  np«-l 

pr 1 nt * » • br anc h  no* from* to* fault  nolif  s.c  is  to  be  simulated)' 
read**br(k )fnfrom( k ) «ntolk  )  *nft (no) 
e l sei f (ch (k ) ( 1 : 1 ) , eq.  i  j >  then 
n  J  =  n  j* 1 
np=no*l 

pr i nt * » ' br anc h  no  * f r om * t o « va lue  of  current  source' 

read**br( k) «nfrom( k) *nto( k) *  cur ( n  j ) 

else 

go  to  100 
endl  f 
continue 

pr  int * »"ambig  set  range  t" 
reao* » amr 

print*»'file  name  where  data  is  to  be  stored:* 
readl 5* * ( a?0 ) * >  filnam 
open(unit=7»file=filnam*status=,new') 
rewind  7 

wr i te ( 7* ' ( "the  first  two  lines  are  for  program  use*  DO  NOT  dele 
te  them")*) 

write(7«»(14i3*el0.3)»)  ne*nl«nnc*nd»nr*nq»ncc*ncv»nvc«nvv*nvt* 
nno»nc  *n j  *amr 
write(7,'(a80) * )  title 

write(7*»("  branch  from  to  control  value  fault  for 
program") ' ) 

write(7*'("type  no  node  node  oranch  no.  use 

only") * ) 
np  =  0 
ne  =  0 

do  120  i=l«k 

i f ( ch ( i ) ( 1 : 1 ) • eo . f e )  then 
ne  =ne ♦ 1 
np=np* 1 


-  30  - 


write (7*250)  ch(i>*br(1)«nfrom(1)*nto(1)«oat(ne)*1e 
endl  f 

120  continue 
nt  =  0 

do  130  1  =  1  ♦  k 

if (ch(1)(i:i). ea.il)  then 
nl^nl+l 
np  =np ♦ 1 

wr1te(7»260)  ch(i)*br(1)*nfrom(1)»nto(1)*nft(np)*1e 
end  i  f 

130  continue 
nnc  =  0 

do  140  i  =  l«k 

i f ( ch< 1 ) ( 1 : 2 ) • eq • i nc  )  then 
nnc  =nnc ♦  1 
np: np> 1 

write( 7*260)  ch(1)«br(1)*nfrom(i )*nto(i)*nft(np)*ie 
end  i  f 

140  continue 
nd  =  0 

do  150  1  =  1  *  k 

i f ( ch ( i ) ( 1 : 1 ) • ea. 1 d )  then 
nd=nd* 1 
np=np*l 

write( 7*270)  ch<1)»br(i)*nfrom(i)*nto(i)*1e 
end i  f 

150  continue 
nr  =  0 

do  160  i  =  1 * k 

iftch(i)(i:i).eq.ir)  then 
nr =nr*  1 

write<7»250)  ch(i>*br(1),nfrom(1)*nto(1)*x(l»nr>*1r 
en di  f 

160  continue 
nq  =  0 

do  170  i  =  1  *  k 

i f ( chi i ) C 1 : 1 ) .ne* i g)  go  to  170 
nq  =  ng-*- 1 

write(7*250)  ch<i)*br<1)»nfrom(1)*nto(1)»x(2»ng)»ig 
170  continue 
n  c  c  =  0 

do  180  i  =  1  *  k 

i  f  (  ch  ( i >  <  1 : 2 )  «ne •  i  c c  )  go  •■o  180 
nc  c^nc  c* 1 

write( 7*280)  ch(i)*br(i)*nfrom(1)»nto(1>»icoi«t(l*ncc)*x(',*ncc)*1cc 
180  continue 
n  c  v  =  0 

do  190  1  =  1  *  k 

if <ch(i )<i:2).ne.1cw)  go  to  190 
ncv=ncv*l 

write<7*280)  ch(i ) »br<1)*nfrom<1 ) »n to  ( 1 )»1cont(2*ncv)*x<4 »ncv)i lev 
190  continue 
nv  c  =  0 

do  200  1  =  1 »  k 

i  f  (  ch<  1 )  (  1 :  2)  .ne.  i  vc)  qo  to  '500 
nvc:nvc»l 

wr1te<7»280)  ch(1)»br(i)»nfrom(1)»nto(1)*1cont(3*nvc)*x{c.  «nvc  ) «  ivc 
200  continue? 
n  vv  =  0 

do  205  i  =  1  »  k 


1 f C ch ( 1 ) ( 1 : 2 ) < ne • 1 v v )  go  to  205 
nwv=nvv»l 

wr1te<7,280>  ch<1>»br(1),nfrom<1).nto<1)»1cont('  nvvWxt'jtnvvWivv 
205  continue 
n  v  t  =  0 

do  210  1  =  1. k 

1  f  (  ch  C 1 )  (  1 : 2  )  »r>e.  1  v t )  go  to  210 

nvt=nvt*l 

np=np*l 

wr1te<7,270)  ch(1)»br( 1).nfromi1)»nto(1).11 
210  continue 
nno  =  0 

do  215  1  =  1.  k 

1 f ( ch ( 1 ) ( 1 12) »ne» Ino )  go  to  215 

nno=nno*l 

np=np* 1 

wr1te(7.260>  ch(1),br<1),nfrom(1>  »nto(1).nft(np)»11 
215  continue 
nc  =  0 

do  220  1  =  1, k 

If ( ch( 1) ( i: 1) .ne.lc)  qo  to  220 

1 f ( ch ( 1 ) ( 1 : 2 ) .eo. 1 c c .or.ch ( 1 ) ( l: 2 ) «eq . 1 cv )  go  to  220 

nc  =  nc*  1 

np=np»l 

wr1te(7»2fc0)  ch(1).br(1).nfrom(1).nto(1).nft(np).11 
220  continue 
n  J  =  0 

do  225  1  =  1,  k 

1 f (ch( 1 ) ( 1 :i) .ne.1 J  )  go  to  225 
n J=n 1 
np=np ♦ 1 

wr1te(7,250)  ch(1>,br(1),nfrofflC 1 ) «  nto( 1),cur(n])«11 
225  continue 

250  formate a 4,14.215, 7x,el0«3,9x«a2> 

260  format ea4, I4,215»17x,13,fcx,a?) 

270  format<a4,l4,215.26x,a2) 

280  formate  a 4, 14,315, 2x,el0«3,9x,a2) 

stop 

end  i 


progam  :  hafdlc 


hafdlc  Is  a  program  for  simulating  hard  faults  in  piecewise 
linear  analog  circuits  and  generating  a  fault  dictionary, 
fault  simualtlon  Is  accomplished  by  a  complemntary  pivot  alo- 
orlthm  for  solving  a  subset  of  tne  circuit  equilibrium  eouat- 
1ons«  which  are  formulated  only  once  In  the  beginning  of  th» 
program  as  detailed  In  vol.l  of  "fault  diagnosis  of  nonlinear 
analog  circuits".  a  dictionary  Is  finally  contained  In  two 
tables  of  testnode  voltage  ranges  and  numeric  fault  codes. 

system:  cdc-6600*  dual  macet  mnf  compiler. 

programmer!  yassin  elcherif 


Input  data: 


a  special  format  file  "data"  must  be  available  In  the  pfiles 
storage  of  the  same  user  Id  specified  In  the  batch  job  card, 
see  program  "cctl"  for  preparing  the  data  file. 

Internal  variables 


bat 

vt 

cur 
nf  t 

header 
br 
type 
nf  rom 

nto 

value 
1  cont 
ch 

kslm 
1  s  1  m 


vector  containing  battery  values 

vector  containing  test  node  voltages  for  the  currently 
simulated  fault. 

vector  containing  values  of  current  sources, 
vector  containing  fault  numbers 
vector  containing  the  problem  title 
vector  containing  branch  numpers 
vector  containing  branch  types 

vector  containing  the  source  node  of  the  cor responding 
branch 

vector  containing  the  destiny  node  of  the  corresoonalng 
branch 

vector  containing  element  values 

vector  containing  control  branches  of  dependent  sources 
array  containing  branch  names  ( 4  characters  each) 
vector  containing  port  branches  to  be  simulated  as  faulty 
under  the  same  fault  number  (multiple  fault  cases), 
vector  flag  for  to  Identify  faults  that  have  been  simulated 


l Imltat Ions 
nan  number  of  nodes 
max  number  of  brances 
max  number  of  batteries 
max  number  of  Independent 
max  number  of  faults 
max  number  of  diodes 
max  number  of  test  nodes 
max  number  of  faults  to  b 


40 

90 
3  0 

cuurent  sources  30 
30 
•50 
30 

simultaneously  simulated  10 


error  messages 

there  are  two  cases  when  a  fault  cannot  be  simulated  • 


in  the  first  case  a  zero  pivot  Is  encountered  In  reducing  tfi-» 
equations  to  formulate  the  complementary  problem.  the  fiult 
will  oe  skipped  ana  the  followinq  error  messeqe  will  be  orinted: 
"fault  no:C  x  3 

zero  pivot-nonzero  c o  l ,  comp l em en t ar y  problem  cannot  be  formulated 
in  the  second  case  the  lemke  complementary  pivot  algorithm  will 
terminate  In  a  ray.  the  fautt  will  be  skipped  and  the  followinq 
messege  will  be  printed: 

"fault  no:Cx] 

no  complementary  solution  Iteration  notxD" 
debugging 

to  print  out  the  complementary  problem  and  the  diode  currents 
and  volta?eS  for  every  fault  simulation  set  the  variable 
1 pdebug=  n  line  2£>0 

to  print  the  same  Information  for  the  nominal  case  only 
set  1pfl=l  In  line  195 


program  hafd1c<data,tape5=data»outout»tape6=output) 
dimension  bat(30)*vt<30>»nft(9Q>»cur(30> 

common/00it/am,q,w,z,ll*b,nll,nl2»a»nel,ne2»1r»mbas1s»1zr 
dimension  am(50«50)»b(50»50)»q(50),a(50),mbas1s(100) 
dimension  w(50)«z<50) 
dimension  nff(30> 

Integer  aa , header ,or « t ype 

common  aa<40«90),ans<90,180>»header(320>»nfrom(90>»nto(90> 
common  br(90)«type(90),vdlue(90),1cont(90) 

Integer  cn<2»90),11,jj»1s1m(90>,lcs1m(10> 

do  1  1=1,90 
nf t ( 1 ) =0 
rewind  5 

read(5»101)  (  header ( 1 ), 1  =  1 , 8 0 > 

read(5,102)  nb,nl,nnc,nd,nr,ng,ncc,ncv,nvc,nvv,nvt,nno,nc, 
read(5 , 101 )  ( header ( 1 ), 1 = 1 , 80 ) 


read<5,101)  ( header ( 1 ), 1 = 1 , 8 0 > 

read(5»102)  nb«nl,nnc,nd,nr*ng,ncc,ncv,nvc*nvv,nvt«nno,nc,nj*amr 

read<5 , 101 )  ( header ( 1 ), 1 = 1 , 80 ) 

wr1te<6,1011)  <header(1),1=l,80) 

wr1te(6,107) 

wr 1 te<  6, 108  ) 

read(5,101)  ( header ( 1 ), 1 =1 , 80 ) 
read(5,101)  < head er ( 1 > , 1 =1 ♦ 80 ) 
np  =  0 
1  =  1 

while (l.le.nb)  do 
np=np*l 

read(5»103>  ch ( 1 , 1 ) , ch<  2, 1 ) , br ( 1 ) , nf rom( 1 ) *n to ( 1 ) ,bat ( 1 ) ,  t  yoe ( 1 ) 
wr1te(6«1031)  ch(l,1),ch(2,1),br(1),nfrom(1),nto(1),bat(1),tyoe(1) 
U1*l 
endwh lie 
n  =  nn*n  l 

whlle(l.le.n)  do 
no-nD» 1 

read<5,104)  ch(l,1),ch(2,1)»br(1),nfrom(1),nto(1),nft(np),type(i) 
wr1te(6,1091)  ch(l*1),ch(2,1),br(1)»nfrom(1),nto(1),nft(np), 

♦  t  ype<  1 ) 

1  =  1*1 
endwhlle 
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n=n*nnc 

wh 11 e (1 • le. n)  do 
np=np*l 

read<5»104)  ch(l*1)tch(2t1>*br(1)«nfrom(1)«nto(1)«nft(no)«tyoe(i) 
wr1te(6*1041>  ch(l»1).ch(2»1)*br<1)*nfrom(1)into(1)*nft(np)» 

♦type ( 1 ) 

1  =  1  +  1 
endwh  1  le 
n=n+nd 

wh 1 le( 1 • le.n )  do 
np  =  np+l 

read(5»105)  ch<l»1)*ch<2»1}»br(1)tnfrom(1)»nto(1)»tyoe(i) 
wr1te<6*1051)  ch(l*1)ych(2»1)»br(1)tnfrom<1)»nto<1)»type(i) 

1  =  1  ♦  1 

endwh 1 Le 

n=n+nr+nq 

wh 1 le< 1 • le .n )  do 

read(5«10  3)  ch(l«1)tch(2t1)*br(  1)»nfrom( 1 ) *  n  t  o  < 1)«value( 1 ) » t  y  o  e  ( 1) 
wr1te<6»1031)  ch(l*1)«ch(2i1)«br(1)infron>(1)tnto(1)«value(l  )  * 
♦type( 1 ) 

1  =  1  +  1 
endwhl Le 

n=n+ncc+ncv+nvc+nvw 
whl Le( 1 . le.n  )  do 

read(5*10£)  ch(l»1)«ch(2»1)«br(1)tnfrom(1)»nto(1)»1cont(1)» 

♦  vaLue( 1 ) *type(  1 ) 

wr1te(6*1061>  ch  ( 1 » 1  )»ch(2*1)ybr(1)ynfrom<1)ynto(1)*1cont(i)« 
♦valuet 1 ) «  t ype(  1 ) 

1=1  +  1 
endwhl le 
n=n*nvt 

whl le( 1 . le. n)  do 
np=np+l 

read(5*105)  ch<l»1>ych<2*1>ibr(1>tnfrom(1)tnto(1)«type<1) 
wr1te<6»1051>  ch(ly1)tch(2y1)fbr(1)»nfrom(1)ynto(1)*type(1) 

1  =  1  +  1 

endwhl le 

n:n*nno+nc 

whl let  1 . le.n  )  oo 

np=np+l 

read(5tl04)  ch  < 1  *  1 > » ch ( 2» 1 >  *br ( 1 >  »nf rom < 1 ) »nto( 1 >  *nf t (no ) ♦ t  ype ( i  ) 
wr1te(6»1041)  ch.  <l«1>«ch(2y1>»br<1)»nfrofr<1)»nto(1>inft<np)» 

♦type( 1 ) 

1  =  1  +  1 
endwhl le 
n=n*n J 
k  =  0 

whl le( 1 . le.n)  do 

np=np+l 

k=k  +  l 

read(5tl03)  ch(l»1)*ch<2»1)»br(1)tnfrom(1)*nto(1>»cur(k>*tyoe(1) 
wr1te(6»1031)  ch(l*1)«ch(2«1)«br(1)«nfron(1)«ntc(1)«cur(k)«tyoe(1> 
1  =  1  +  1 
endwh lie 

101  format(80al ) 

1011  format(lhl»30x»80al> 

102  format(1413*el0.3> 

103  format(a2ta2*14t2l5*7x«el0«3'9xya2> 

1031  for(nat(2x»a2»a2yl4»2l5t7xfel0.3»9x»a2) 

104  format(a2ta2»14»215*17xt13»6x»a2> 


1041  format (2*  ta2fa2»l4,2l5,17x,13»6x,a2) 

105  f orma t ( a2 ♦ a2 ♦ 14 »2  15 1 26 x »  a2 ) 

1051  format <2x»a2»a2« 14*215 . 26x,a2) 

106  format(a2»a2«14»315,2x»el0.3*9x»a2> 

1061  format(2x»a2*a2«14«315»2x*el0.3«9xta2) 

107  format(*  branch  from  to  control  value  fault’) 

108  format(’name  no  node  node  branch  no’) 

c 

ndv=nd*nvt 

nl=nb*nl*nnc 

ntb  =  np*nr*ng«-ncc*ncv*nvc*nvv 
c  obtain  hybrid  matrix 
c 

call  hybr Id ( 1 1 ♦ J J ♦ 0 m t b ) 

Ipf 1=1 
c 
c 

c  obtain  number  of  non-port  branches 
c 

11=11-1 

n=n-i 

c 

c  form  and  solve  the  complementary  problem  of 
c  the  nominal  case 
c 

do  4  1=l«ndv 
q< 1 ) =0.0 
do  2  J=l»nd 

2  am( 1 »  J  >  =  ans ( 1 *nl*1 1 ♦ j  *np*n  1* J  J  ) 
do  3  J  =  1  mb 

3  q(1)=q(1)*bat(j)*ans(1*nl*11»j*np*J}> 
n=nl*nd+nr»ng+ncc+nc  v*nvc*nw*nvt  *nno*nc 
If(nJ.ne.O)  then 

do  303  j  =  1 m j 

303  q( l)=q(1)*cur(j)*ans( 11+nl+ltJJ^np^n^J) 
end  i  f 

4  continue 

1 f (ipf  l  .eq.  1 )  go  to  ?01 
wr1te(6» 19) 
ndg  =  0 
1=1 

2012  ndg  =  m1n0  (nd»ndg*9 ) 
do  2020  1  =  1  md 
1  f ( ndg. 1 1 .nd )  then 
wr1te(6»20>  ( am( 1 « j ) « j =1 »ndq) 
else 

wr1te(6*20)  (am(1«J)»J=l*ndg)*o(1) 
end  1  f 

2020  continue 
1  =  1*9 

1 f (ndg. I t .nd )  no  to  2012 
201  call  lemke(nd) 
c  zero  out  vt(1) 

c  calculate  and  print  out  node  voltages  of  the  nominal  case 
c 

do  5  1=l*nvt 

5  vt(1)=0.0 

do  7  1=l»nvt 
vt(1)=vt(1)-q(1*nd) 
do  7  J=l»nd 


7  vt ( 1 > =vt ( 1 > -am ( 1 ♦nd* J ) *z  ( J  )  ... 

do  161=1 *nwt 
16  ans( 1  *  1 >  =vt (  1 ) 
n3=nl 

n4=nl+nd^nr+ng*ncc*ncv+nvc+nvv 
wr 1 te (6 • 10 ) 

wr1te(6*8)  (ch(l*1en3)*ch(2»1^n3)»w(1)*z(1)»1=l*nd> 
wr1te( 6* 11 ) 

wr1te(6,9)  (ch(li  1*n4) ieh(2»1  +  n*),vt ( 1) ii:l»nvt ) 

10  format (//10x * ‘diode  cur r ent ' *  1  Ox » • d 1  ode  voltaqe‘> 
c  format  (/5x*2a2,el0.3»lr)x,el0. 3) 

11  f  ormat  ( //2 x  * nom  1  na  L  test  noce  wo 1 1  ages  ’ //5x  * ‘node  *  *  r  x  * •  vo 1 1  age  '  > 

9  f ormat ( / 5x ♦ 2a 2* 4 x , e 1 0 . 3 ) 

19  f ormat ( /2 x» como lement ary  problem* //5 x *• coef f  matrix*  last  col-i1 > 

20  format(/10(lx»elQ.3)> 

18  f ormat ( //2x , ‘no  of  Iterations  =  ‘»13> 
nbl=nb+l 
nnf  =0 
zero=l.e-8 
1 debug=0 
do  400  1=1*90 
400  1s1m(1)=0 

nnf  =  1 
nf  f ( 1 ) =1 
c 

c  loop  to  simulate  faults 
c 

do  90  nf=nbl»np 
If(nf.ge.n)  go  to  90 
1 f ( n f t ( n f ) • eq. 0 )  go  to  90 
1 f ( 1 s 1m (nf ) .eg. 1 )  go  to  90 
do  320  1=1*10 
320  kslm(1>=0 

m  =  0 

do  330  k=nf«np 

1 f < nf t C k ) .ne.nf t (nf ) )  go  to  330 
m  =  m*l 
ks1m( m) =k 
1 s1m( k ) = 1 
330  continue 
1 f lag=0 

1  f (nf «gt «nl .and.nf . le • ( nl  +  ndv ) >  go  to  90 
c 

c  simulate  only  faults  Identified  In  the  Input 
c 

150  do  23  1=1, ndv 
0(1 >=0.0 
do  21  J=l*nd 

21  am(1*J)=ans(1mnlm11»Jmnp*nleJJ) 
do  22  J=l*nb 

22  q( 1 >  =  q( 1 > ♦bat ( J  >  *ans ( 1 enl* 1 1 »  j*np* J  J ) 
n=nl*nd*nr+ng*ncc*ncv*nvc+nvv»nvt*rmo*nc 
1 f (n J »ne  .0  >  then 

do  322  J  =  1  *  n j 

322  o(1)=q(1>»cur(J)*ans(11^nlm1*JJ*nomn»J> 
endl  f 

am( 1 ,nde  1 >  =a ( 1  ) 
do  340  l=l«m 
k  =  ks1m(  l) 

340  a*(1*nd»lml)=ans(11»nlm1*JJ>np*k) 


23  continue 
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c 

c 

ndvm=ndv+m 
ndwl  =  ndv<,l 
do  345  i=ndvl»ndvm 

345  am< 1 » nd* 1 ) =0 . 0 
c 

do  360  i=ndvl»ndvm 
k  f=ksim< 1-ndv) 
do  346  j=l*nb 

346  am<1*nd+l>=am(i*nd«-l)*bat<j)*ans<ii«-kf*jj«-np  +  j> 
do  347  j=l*nj 

347  am<i»nd-*-l>=am(i»nd-*'l>+cur<j>*ans(11  +  kf»Jj  ♦no+n+ j  ) 
c 

c 

do  348  j=l*nd 

348  am(i»J)  =  ans(11^kf*JJ'*'np*nl*i) 
c 

c 

do  349  1  =  1 »  m 
k=ks1m(l) 

349  am( iind  +  l+l )=ans(ii*kf » j j  »np»k) 

360  continue 

do  390  l=l*m 
pivot=am<ndv  +  l»nd-»'l  +  L) 

1 f ( abs(p1  wot )  .  It. zero )  then 
do  370  i=l»ndvm 

if(abs{am(i»nd«-l*l>>.qt.zero>  then 
print  361 *nf t (nf  ) 

361  f ormat ( // • f aul  t  no : • « 1 2* *  zero  pivot*  nonzero  column  in  the** 

♦  •hybrid  mat r i x • / *  comp l emen t ar y  problem  cannot  be  formulated*) 

write (6, 362)  ( (am ( i k » j k ) , ) k  =  1 *6 ) ♦ i k  =  l . 6 ) 

362  f ormat < 6 ( 2 x *  el  0 • 3 ) ) 
go  to  90 

end  1  f 

370  continue 

if(m.eq.l)  then 
nnf  =  nnf«-l 
do  36  i=l*nvt 
36  ans (  i  *nn f )  =  ans ( i » 1 ) 
go  to  90 
end  1  f 
go  to  390 
endi  f 
ndm=nd«-m 
do  372  j=l»ndm 

372  am(ndv^l*J)=am(ndv*-l»j)/pivot 
do  374  i=l»ndvm 
do  374  J=l*ndm 
1 f (  i  • eq • (nOv+ 1 ) )  go  to  374 
am(i*j)=am(i»j)-am(ndv*l*J)*am(i»nd>l'»l) 

374  continue 
390  continue 
c 

do  28  i=l*ndv 
28  q ( i  )  =am < i  ♦nd  +  1  ) 

i  f (  1  f  lag .eq « 0 )  go  to  209 
1 f ( 1  debug  .eq.  0  )  go  to  209 
print  210 


210  f orma t ( 2* • • bef ore  attemptinq  to  solve  by  lemke  alqorlth^ 
♦  the  prob lem  was  » ) 

go  to  211 

209  call  lemke(nd) 

211  If (1r.ne.1001)  go  to  206 

1 f< 1  debug. eq. 0 )  go  to  206 
If ( 1 f lag.eq.l)  go  to  3011 
wr1te(6»15)  nft(nf)»nf 
wr It e( 6. 19 ) 

3011  1=1 
ndg  =  0 

3012  ndg=m1nQ(nd.ndg+9 ) 
do  3020  1=l.nd 
If(l.eq.l)  then 

wrl te( 6* 301 3 )  (  c h ( 1 . n Z* } ) . ch ( 2 « n3«- ] ) . j  =  1  * nd ) 
endl  f 

3013  format(10(lx.2a2»4x)  ) 

1  f ( ndg. 1 1 .nd )  then 
wr1te(6»20)  (  amt  1 ♦ J ) » j  =  l »ndg ) 
else 

wr1te(6.20)  (am(1.j).j=l«ndq).q(1) 
endl  f 

3020  continue 
1=1*9 

1 f (ndg. I t .nd)  go  to  3012 
wr1te(6.18)  11 
wr1te(6«10> 

15  format (//2x» ' fault  case  no!  *» 12»5x *• column  no!  '.12) 

1 f ( 1 f lag .eq. 1 )  go  to  206 
1f( ldebug.eq.0)  go  to  206 

wr 1te(6»8)  ( ch  ( 1 . 1  *n  3 >  »ch  (  2. 1  *n 3 )  .w  ( 1  > .  z  ( 1 )  ♦  1  =  1  .nd ) 

If lag=l 
go  to  150 
c 

c  calculate  node  voltages  under  fault  condition  nft(nf) 
c 

206  do  35  1=1. nvt 

35  vt(1)=0.0 

do  37  1=l.nvt 
vt ( 1 ) =- o( 1*nd) 
do  37  ]=l.nd 

37  vt ( 1 ) =vt ( 1 ) -amt 1*nd» J ) *z ( J ) 

1 f ( Ipf l .eq. 1 ) go  to  207 

w  r 1 te ( 6 «  39 ) 

38  format(//5x»' node*  * 5x .  • voltage*) 

wr1te(6.9)  (ch(l*1*n4)»ch(2»1«-n4)»vt(1)»1  =  l»nvt) 
c 

c  store  the  node  voltages  of  the  current  fault  case  In  the 
c  f-vt  table 
c 

207  1  f (  1r.eq.1001)  go  to  90 
nnf  =nnf ♦ 1 

nf f (nnf ) =nf t (nf  > 
do  39  1=1. nvt 

39  ans( 1»nnf)=vt( 1) 

90  continue 

c 

c  print  out  f-vt  table 


nfmax=nnf 


wrl te(6»41 ) 

41  f ormat ( //2x » • f au 1 1 -v t  table*) 

1  =  1 

110  =  0 

40  1 10  =  am1n 0< 1 10*1 0  inf ma x  ) 

wr1te(6*42)  ( nf f ( j ) * j  =  l « 110 > 

42  f ormat ( //2x » *  fault  no » « 2 x «  1 0 (1 2 *9 x  ) ) 
do  43  1=l»nvt 

43  write(6#44>  cM  1  ♦  1  *n  4  )  «  ch  (  2 »  1  *n 4 )  t  ( ans  ( 1 1 J  >  ♦  j  =1 1 1 1  0  ) 

44  format(//2x»2a2t4x»10(el0«3»lx)  ) 

1 f ( 110 • ge* nf max)  go  to  111 
1=1*10 

go  to  40 

111  call  ambset (nvt «nf ma x i ch * n4 , amr ♦ nf f ) 
end 


c 

subroutine  hybrid  for  obtaining  the  hybrid  equations  of  the 
c  n  port  network, 

e 

c  Input  and  common  variable0, 

c 

c  a 

c  ans 

c 

c  nportl 

c  1 1 

c  debug 

c  ntb 

c 

c- - 

c 

subroutine  hybric<nportl»11«1ornt»ntOtdebuq) 

Integer  a»nfrom(90)»nto(90)»1cont(90)»type(90)»dcol(90) 

Integer  1count(?)»header(32Q)»br(90)*rbr{90>»ansrow»ansrol 

Integer  cv«etr»cc»vv»vc»from,to»c«coun»oeg1n»temp»st»tn 

Integer  tp * por t t oebua« 1 s t p « q» v 

real  value(93)«f3(90»40>«f6<90t40>»dnsC90tlBO> 

common  a ( 40  * 9 0 )  ians»header»nfrom*nto 

common  br « t ype • va lue» 1 cont 

data  e  t1s«rtg»vv»cvtccivctsttc»v/2h  e»2h  1»2h  r»2h  q» 
12hvv*2hcv»2hcc»2hvc»2hst*lh1»lhv/ 
data  Isto/lhs/ 
c 

c  max  circuit  conf 1 gurat 1  on  Is  40  nodes  and  80  branches 
c  max  of  40  elements  In  any  category  (tp»tn*ln*lp) 
c 

1  nbr  =  0 
nnodesQ 

c 

c 

c  zero  out  a  matrix 
do  2  1 =1 ♦ 40 
do  2  J=lt90 

2  a ( 1  *  j ) =0 

c  read  In  data  and  fill  a  matrix 
debug=0 
do  3  k=l»ntb 

c  store  entries  Into  a  matrix 
f  rom=nf rom( k ) 
to=nto(k) 

1f(from.re.0)aCfrom«k)=l 
1f(to.ne.0)a(to*k)=-l 
nbr =max  0 ( nb  r  «  k  ) 

3  nnode=mdxO(nnode«nfrom(k)«nto(k)) 
zero=l.e-fl 

1 f ( debug.ne. 1 )  qo  to  5 
c  print  the  a  matrix  for  debug  run 
wr 1te(6» 505) 
do  6  1=l»nnode 

6  wr1te(6t506)(a( It J)» J=l*nbr) 

5  do  7  1=l»nbr 

dcol (1  )  =  0 

7  rbr< 1 ) =0 
c 


array  containing  Incedence  matrix 

array  containing  the  hybrid  matrix,  In  the  reduced 
echelon  form 

starting  row  of  the  hybrid  matrix  In  ans 
starting  column  of  the  hybrid  matrix  In  ans 
flag  for  printing  debugging  outout 
total  number  of  branches  In  the  network 


O'  eo 
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c  determine  elements  maklnq  up  the  tree 
c 

call  f t reel nnode* rbr ♦ dco l ) 
c 

c  reorder  a  matrix  Into  four  classes 
c 

c  1«  tree  port  branches(tp) 

c  2«  tree  non-port  branches(tn) 

c  3«  link  non-port  branches(ln) 

c  4.  link  port  branches(lp) 

c 

c  dcol  contains  ordering  of  a  with  tree  branches  In  leftmost  columns 
J  J=nnode*l 
n  =  l 

do  8  j=l«nnode 
m=dcol ( ] ) 
do  9  k  =  n»m 
1  f ( m. eq. k ) go  to  8 
dcoKjJ  )=k 

JJ  =  JW 

contl nue 
n  =  m>  1 

1f<  JJ  •eq»ncol)oo  to  10 
do  11  1=n»npr 

11  dco l  <  1 )  =  1 

c  reorder  dcol  Into  four  classes 

c  lcount(l)  marks  last  port  column  of  tree  branches 
c  1count<2>  marks  last  non-port  column  of  link  branches 
10  1  count ( 1 )  =  1 

1 t2=nnode 
1=1 

13  do  12  m  =  1 1 1 1  2 
mm=<nbr«-l)*(1-l)M<3-<2*1>>*m) 

1  tem  =  dcol( mm) 

1 f ( t /pe< 1 1 em ) .ne.e ,and • t ype ( 1  tern ) »ne. 1 s )go  to  12 

1  tern  1  =  1  count ( 1  ) 

dco 1 1  mm ) =  dco l ( Iteml ) 

dco l  ( 1 1 eml ) =1  tern 

1 c  ount ( 1)  =  1count ( 1 > ♦ 1  —  < < 1-1 ) *  2 ) 

12  continue 
1f(1.eq.2>  go  to  14 
1count(l)=1countCl)-l 
1  count ( 2  >  =  nbr 

1  t2  =  nbr-nnode 
1  =2 

go  to  13 

c  reorder  the  a  matrix  and  the  original  label  vector  to  correspond  to 
c  the  reordered  dcol 

14  nn- 2 
n=  1 

beg1n=l 

coun=0 

15  1 tem  =  dcol (n  ) 

1 f< 1 tem.eq. begin) go  to  18 

1 t emp=or ( n) 

br<  n) =br ( 1 1  em) 

br ( 1 1 em)  =  1 1 emp 

do  17  J=l«nnode 

t  emp  =a< J in ) 

a< J»n)=a< J» Item) 
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a  ( j « 1 1  em )  =  t  emo 
coun=coun*l 
dcol(n)=-dcot(n) 
n=1 tem 
go  to  15 
16  dcol (n )  =  -dcot (n  ) 

1 f < coun.eq. <nbr-l )  )qo  to  18 
do  19  1=nn»nbr 
1 tem  =  dco  l  ( 1 ) 

1 f ( 1 t em »eq • 1 ) go  to  20 
1 f ( 1 tem. I t • 0 ) go  to  19 
beg1n=1 
n  =  1 

go  to  15 

20  coun=coun+l 

dcol(1>=-dcot(1) 
nn  =  1 

19  continue 

18  do  22  n=l»nbr 

22  dcol (n ) = 1 abs ( dco l ( n ) > 
c 

c  reduce  reordered  a  matrix  to  row  echelon  form 
c 

call  laech ( nnode «nbr ) 
c 

c  back  substitute  a  matrix 
c 

do  23  1=2»nnode 
lrow=1-l 
do  23  J=ltlrow 
1  f col  =  1 

1  temp=a( J» Ifcol) 
do  23  k  =  1  »nbr 

23  a ( J t k ) -a ( J « k ) -a ( 1 « k ) *  1 1 emp 
c 

c  formulate  the  element  characteristics 
c 

c  tp  Is  number  of  columns  In  fl  and  f*5 

c  tn  Is  number  of  columns  In  f2  and  f& 

c  In  Is  number  of  columns  In  f3  and  f7 

c  Ip  Is  number  of  columns  In  f4  and  f8 

tp=1count( 1 ) 
tn=nnode-1 count (1) 
ln=1count(2)-nnode 
lp=nbr-1count (2) 
port=tp>lp 
nport=tn* In 
ansrow=nbr 
anscol  =nbr*por t 
wr1te(6*507  > 

If(tp.eo.O)  go  to  24 
wr1te(6»508)  ( br ( 1 ) ♦ 1  =  i , tp  ) 

24  J=tp*l 

1 f ( tn.eq.O )  go  to  25 
wr1te<6»509 )  (br ( 1  )* 1  =  J , nnode > 

25  J=nnode*l 
J  J  =nnode*ln 
If(ln.eq.O)  go  to  28 
wr1te(6«5l0 )  ( br ( 1 > *  1 = j ,  1 J > 

J  =  JJ*1 
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1 f ( lp.eq.0 )  go  to  28 
ur1te(6«5ll)  ttrU  >.  i  =  1  fnbr> 
c  zero  ans  matrix 

28  do  27  1=liansroM 
ao  27  j=l*anscol 

27  ansi t*J)=0.0 

do  29  1=l«nport 
do  30  J  =  1  *  t  n 

30  f6(1tj>=0 

do  29  5  =  1*  In 

29  f 3( 1 » J  )  =  G 
kount  =  1  count ( 1 ) 
k  =  0 

!=1 

do  31  1=l«nbr 
1tem=br<  1) 

31  rbr  C 1 tem )  =  1 

1 f C debug. ne. 1 )  go  to  32 
wr1te(6*512>  tpttn«ln«lp 
wr1te(6»513)(br(1)»t=l»nbr) 

32  kount=kount*l 
mm=dcol( kount ) 

1 temp=1cont (mm) 

1temp=rbr(1temp) 

1 tl=pcrt»  J 

1f(type(fflm).eq.g.O''.type(mm>.eq.vc»or.type(mm).eq.cc) 
lgo  to  33 

c  voltage  source  type 

If  ( kount .gt .nnode >go  to  34 

c  f  2 

1t2sln*1 

ans (1tl*1t2)=l. 

If (tyye(mm) .eq.cv)go  to  35 
1 f ( t y pet mn) . eq. v v) go  to  36 
f6(J»J)=-vaiue(mm) 
go  to  37 
34  k  =  k*l 

f 3 ( J*k)=l. 

1 f ( t ype < mm ) «eq . c v > go  to  35 
1 f ( type( mm) .eq.vv) qo  to  36 

c  f  7 

ans( 1tl«k)=-value(mm) 
go  to  37 

c  current  source  type 

33  1 f (kount .gt  .nnode)go  to  38 
f6< J» 5  >  =  1 • 

1 f < t ype (mm ) .eq  •  vc  )qo  to  36 
1 f ( typet mm) .eq. cc)qo  to  35 

c  f  2 

1t2=ln>J 

ans(1tl*1t2)=-value(mm) 
go  to  37 
38  k  =  k*l 

c  f  7 

ans(  1tl»k)=l. 

1 f ( type(mm) .eq.vc )qo  to  36 
1 f ( type ( mm) .eq. cc > go  to  35 
f3(J»k)=-value(mm) 

37  J=Jel 

1  f( kount .ne . 1  count (?)) go  to  32 


go  to  39 

c  current  controlled 

35  1 f < 1  temp «gt • t p )go  to  40 
c  f5 

1t2=nport*1temp 
ans(1tlt1t2)=-value(mm) 
go  to  37 

40  If < 1temp.gt.nnode)go  to  41 

I I  =  1 1  emp-tp 
f6(J*1t)=-value(mm> 
go  to  37 

41  1 f( 1  temp • gt • 1  count ( 2 ) )go  to  42 
1t=1temp-nnode 

c  f  7 

ans  ( 1 1 1 * tt )  =-value(mm) 
go  to  37 

42  1 t = 1 t emp-1 count ( 2) 
c  f8 

1t2=nbr«-tp*1t 
ans(1tl*1t2)=- value (mm) 
go  to  37 

c  voltage  controlled 

36  1 f ( 1 temp.gt • tp> go  to  43 
c  ft 

1 t2=nbr *1  temp 

ans( 1  tl  t 1t2)=-value(mm) 

go  to  37 

43  1 f ( 1  temp .gt .nnode  )  go  to  44 
1 t = 1 temp-tp 

C  f  2 

1t2=ln*1t 

ans< 1tl«1t2)s -value (mm) 
go  to  37 

44  1 f( 1 temp .gt « 1  count ( 2 )  )go  to  45 
1 1=  1  temp-nnode 

f  3  <  J  tit)--*  value  (mm) 
go  to  37 

45  1t=1temp-1count<2> 
c  f  4 

1t2=nport-*-tp*1 1 
ansl1tl*1t2>=- value (mm) 
go  to  37 

39  1 f( debug «eq« 0 )  go  to  46 

If(ln.eg.O)  go  to  47 
c  write  f 3  for  debug  run 
wrt  te<6$514  ) 

I I I  =  1 

48  1 12  =  In 

If C (1t2-1tl)*gt,10)  go  to  49 
1f( 1t2.eq.1tl)  go  to  47 
wrl te  <6 1515 ) 
do  50  1=ltnport 

50  wrl te ( 6 t 516 )  < f 3 ( 1 t J > t J = 1 t 1 t 1 t 2 ) 
go  to  47 

49  1t2  =  1tl«-9 
wr1te<6t515) 
do  51  1=l*nport 

wr 1 te( 6*516 )  <  f  3  <  1  *  )  >  *  J  =  1 1 1 » 1 1 2 ) 
1 1 1 1 2*1 
go  to  48 


51 


47  If(tp.eq.O)  go  to  52 

c  write  f6  for  debug  run 
wrl te  <  6»  517 ) 

1 1 1  =  1 

53  1 t2=tn 

1f<<  H2-1tl).gt.l0>  go  to  54 
If  Clt2.eq.1tl>  go  to  52 
wr1te(6*515) 
do  55  1=l*nport 

55  wrl t e<  6  *516  )  ( f 6 ( 1 . J ) * J  =  1 1 1 ♦ 1 1 2  ) 
go  to  52 

54  1t2=1tl*9 

wr 1 te ( 6  *515  ) 
do  56  1=l*nport 

56  wr1te(6»516)  <  f  6  <  1  .  J  >  * }  =1 1 1  » 1 1  2  ) 

1 tl = 1 t2*  1 

go  to  53 

52  wr1te(6»518) 

call  pr int ( ans co l * ans row ) 
c 

c  zero  out  f6 
c 

46  1 f ( tn.eq. 0 ) go  to  57 

do  58  J=l»tn 
kk  =  tp«- ] 

do  58  1=l*nport 
1 tl=port*1 

If(ln.eq.O)  go  to  59 
c  change  f7 

do  60  k=l«ln 
lk=nnode*k 

60  ans( 1tl» k)=ans( 1tl*k)-< f6< 1* J)*a(kk»lk) ) 

59  If(lp.eq.O)  go  to  58 

c  change  f8 

do  61  k  =  l*  Ip 
lk  =  1  count  <  2 ) ♦ k 

I  t2  =  nbr  +  tp«'k 

61  ans< 1 1 1  *  1 1 2 )=ans(1tl»1t2)-(f6(1»J)*a(kk*lk)) 

58  continue 

c 

c  zero  out  f3 
c 

57  1 f ( In.eq.C )go  to  62 
do  63  J=l»ln 
lk=nnode* J 

do  63  1=l*nport 
1 1 l  =  por t  +  1 

If(tn.eq.O)  go  to  64 
c  change  f2 

do  65  k=l*tn 
kk=tp*k 
1 1 2= In*  k 

65  ans(1tl*1t2>=ansC1tl*1t2)-<f3<1» j > * <-a( kk* lk> > > 

64  1f(tp«eq*0)  go  to  63 

c  change  fl 

do  66  k  =  1  *  t  p 

I I  2=nbr *k 

66  ans<1tl*1t2)  =  ans(1tl*1t2)-<f3<1»  J)M-a<k,lk>  >> 

63  continue 

c 


c  fill  ans  matrix 
c 

62  1  f ( debug .eq • 0 )  go  to  67 

wr1te<6t519> 

call  pr  Int (anscol  tansrow  ) 

67  If  (In.eq.O.or.tp.eq.O >go  to  68 
c  store  dl 

do  69  1=l»tp 
do  69  J=i i in 
k  =nnode* j 

69  ans(1«5)=a(  1«k) 

68  lc=ln*l 
ltempstp+1 

1 f ( 1  temp. gt • por t .or. I c • gt .nport ) go  to  70 
c  store  -d4  transpose 

do  71  1=1temp*port 
JJ=lc*1-1temp«>nnode 
do  71  j  =lc tnport 
1 1  =  J  +  l-  lc*tp 

71  ans<  1 ,  j )  =-a  C 1 1 » j  J ) 

70  If(tp.eq.O)  go  to  72 

c  store  unit  matrix  above  f5 
do  73  1  =  1 1 1  p 
ld=nport>1 

73  ans<1,ld)=l. 

72  If(lp.eq.O)  go  to  79 

c  store  unit  matrix  above  f4 
1 1  =  tp+l 

do  75  1=11 t  por  t 
ld=nport+1 

75  ans(1,ld)=l. 

74  1temp=tp*l 
If =ld*tp 
le= la*l 

1  f  (  1 1  errp  .gt  .por  t  .or  .  le  .gt  •  l  f  )  go  to  76 
c  store  ~d 2  transpose 

do  77  1=1temp.port 
JJ  =  1-1temp«-1count<2)*l 

do  77  J= le»  If 
1 1 =J ♦l-le 

77  ans( It} )  =  -a< 1 1 » J  J ) 

76  le=lf*lp 
l  d  =  l  f  ♦  l 

1 f ( tp. eq. 0. or « l d .gt • l e ) go  to  78 
c  store  d2 

do  79  1=l*tp 
do  79  J  =ldtle 
k=1count  <21 ♦ J-lo 
79  ans( 1« j )  =  a< 1«k  ) 

78  1 f (debug .eq. 0  )  qo  to  90 
wr1te<6.520) 

call  prlnt(ariscoltansroM) 
c 

c  reduce  ans  matrix  to  echelon  form 

call  raech(nbrtanscoltar.'col  •  1 « 1 « ;er o  ) 

1 f < debug. eo. 0 )  go  to  81 
wr1te<6»521  ) 

call  pr Int ( ansco l . ans r om > 
do  82  1 = 1 tnbr 
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do  82  J=l,nport 

1 1  =  nb  r *  1 -  1 

If  (  abs(  ans(  H  *  J  ) )  .  le  .  zero  >  ansC11,J>=0.0 
If  Cans ( 1 1 t J ) .ne.O . ) go  to  83 

82  continue 

83  11=11*1 
c 

c  fill  column  heading  vector  for  final  print  out 
c 

J=0 

If(tp.eq.O)  go  to  84 
do  85  1=l»tp 
1 t=2*1 

headerC 1 t ) =brC  1 ) 
header(1t-l)=c 

12  =  2  * ( port*  1 ) 
headerC 1 2) =br  C 1 ) 

85  header C 12-1  )  =  v 

84  If(lp.eq.O)  go  to  36 
J  =tp 

do  87  1=1, Ip 

J=i*l 

k=  1* 1  count  C  2 ) 

1 1  =  2*  j 

headerC It )  =  brCk ) 
headerC1t-l)=v 
1 2=2*  C  por t*  tp* 1 ) 
header C 1 2)  =  br C  k ) 

87  header  C 1 2-1 )  =  c 

86  1  t  =  4*port 
nport l=nport  *1 
do  88  1=11, nbr 

do  88  }=nportl, anscol 

88  1 f CabsCansC  i ,  J ) ) . le. zero)  ansC1,J)=0.0 
1 fC debug. eq. 0 )  go  to  89 

c  print  final  ans  matrix  for  debug  run 

call  p r nt  1 C 1 1 , npor 1 1 , anscol, 1 1 ,nbr ) 

89  If  C11.eq.nbr)  go  to  90 
c 

c  back  substitute  final  answer  matrix 
c 

1 t 1 =ansrow- 1 1 *1 
1 1 2  =  1 1  *  1 

do  91  1=1t2»ansrow 

c  ans(1rw»1cl)  Is  pivot  element  using  to  zero  elements  above 
1 rw=ansrow* 1 1 2  —  1 
1cl=nport*1 t 1 ♦ 1 t2-1 
1 1 3  =  1 rw-1 

c  J  =  row  zerolnq  out  above  pivot 
do  91  J=11,1t3 
b=ansC  J,  Icl) 

c  k=column  chanqlng  of  Jth  row 
do  91  k=1cl,anscol 

91  ansCJ,k)=ansCJ  ,k)-b*ansC1rw«k) 

90  do  92  1=11, nbr 

do  92  J=nportl, anscol 

92  1 f C abs C ans C 1 ,])). le • zero  )  ansCl,))  =  0.0 
c 

c  print  final  ans  matrix 


1  f < Iprn t »eq .0 )  go  to  499 

call  prnt 1 ( 1t»nportl» an scol»11*nbr> 

return  ■  1 

format  (80al) 

format(lhl/lx»80al//////lx»*network  descr1pt1on*»///lx»» branch 
from  to  element  element  control*/*  number  node  node 

r  type  value  branch*) 

format<313»a2»fl0.3» 13) 

format(2x«13»6x»13»3x»13*7x*a2*4x»el0.3»3xi13) 
f orma 1 1 lh 0// *  *  zero  =  *iel0.3> 
formate///*  a  matrix*) 
f ormat ( IhO  ?  50 1 2  > 
format ( IhO/// > 

f ormat ( 1  ho • t ree  port  branches ' /50 ( 1 x *  1 2 ) ) 
format ( lh0« 'tree  non-port  b r anc he s * /5 0 ( 1 x » 1 2  )  ) 
f ormat < IhO ♦ * l Ink  non-oort  b ra nc he s * / 50 ( 1 x « 1 2 ) ) 
f  o»*mat  <  IhO  *  *  1 1  nk  port  br an ches * /5 0 ( 1 x » 1 2 >  ) 

format < lh0» *tp  =  *»13/*  tn  =  **13/*  In  =  *»13/*  Ip  -  »,13) 
formate lh0»*br*  »50<lx»12) ) 
formate///*  f3  before  zeroing*) 
formate  1 x) 

format(lx»10(ell.4»lx)) 
formate///*  f 6  before  zeroing*) 
formate///*  ans  matrix  before  zeroing*) 
formate///'  ans  matrix  after  zerolno*) 
formate///*  ans  matrix  with  d  values  filled  In*} 
formate///*  ans  matrix  reduced  to  echelon  form*) 


subroutine  ftree(nrow«ncol«1ndcol) 
c 

c  subroutine  ftree  takes  the  matrix  a*  apolles  subroutine  laech  and  fines 
c  the  Independent  columns  In  a  closest  to  the  left,  these  Independent 
c  columns  make  up  the  tree  branches. 


Integer  a « 1 ndco l ( nrow ) « co l « t emp 
common  a(4Q»90) 

1  =  1 

temp=l 

call  1aech(nrou»ncol ) 
c  step  through  rows 
do  1  k=l»nrow 
c  step  throuqh  columns 
do  2  J=temp»ncol 
e  find  Independent  columns 
c  test  If  element  equal  to  one 
If  < a< k ♦ J ) .ne. 1 )  go  to  2 
c  record  Independent  column  number 
Indcoll l>= ) 

1=1*1 
temps  J*i 
go  to  1 

2  continue 

1  continue 

return 
end 
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subroutine  iaec h( nrow » ncol ) 
c 

c  subroutine  iaech  manipulates  matrix  a  into  echelon  form 
c 

c-- — - - 

Integer  a « c ♦ g ♦ gp l us  1 « p * b 
common  a<40»90) 
c  =  1 
9  =  1 

2  do  1  1=g»nrow 

1 f ( a( 1 # c> • eq»0 ) go  to  1 

c  interchange  1  and  g  row  to  get  nonzero  pivot 
1 f ( 1  • eq*g  >  go  t o  3 
do  4  k=c*ncol 
b=a  (1 *k  ) 
a(i*k)-a(g*k) 
a (g  »k  )  =b 
4  continue 

c  normalize  row  to  get  positive  number  for  pivot 

3  1 f < a (g »c  )  .g t  •  0  )  go  to  5 
do  6  k=c»ncol 
a <  g « k ) =-a ( g * k  ) 

1 f ( g« ge • nrow)  return 

c  zero  column  below  pivot 
gplusl=g+l 
do  7  p=gplus 1 tnrow 
b=a ( p t c ) 

1 f ( b. eq . 0 ) go  to  7 
do  8  k=c*ncol 

8  a(p» k ) =-b*a< k ) ♦a(p «k  ) 

7  continue 

9=9*1 
c=c+l 
go  to  2 

1  continue 

1 f ( g. gt .nrow )  return 
c  =c  ♦  1 
go  to  2 
end 
c 

c  - - - - - - - - — - - 

c 

subroutine  raech(m»n»mark»rowl»coll»zero) 
common  1  a (  4 0  *  9 0 ) « a ( 3 0 » 18 0 > 

Integer  c »g »gp l usl » d « rowl*c ol 1 

raech  performs  row  operations  on  a  to  reduce  a  to  echeton  form 

columns  coll  to  mark  are  reduced  to  row  echelon  form  while  the  row 
operations  are  carried  out  on  the  rows  from  mark  ♦  l  to  n 
rows  rowl  to  m  are  reduced  to  row  echelon  form 
g=row  determining  pivot  point  in 
c=column  determining  pivot  point  in 


c  find  the  max  nonzero  element  In  the  c  column  below  and  Includina  pivot 
1  =  0 

tmz  =0 • 0 
do  2  j  =  g*m 

1f(abs(a(J*c)>«le«zero>  a<  J  *  c  >  =  0  •  0 
1 f (abs (a ( J tc  ) )  •  le.tmz  )  go  to  2 
tmz  =  abs( a( J  « c) > 

1=1 

2  continue 

If (tmz.eq.0.0 )  go  to  1 

c  If  the  nonzero  element  is  in  the  pivot  row*  do  not  exchange  rows 
If ( 1 .eq.gigo  to  3 

c  exchange  pivot  row  with  row  having  nonzero  element  In  pivot  column 
do  4  k  =  c*n 
b  =  a ( 1 « k  ) 
a(1*k)=a(q«k) 

4  a ( g*  k ) =b 

c  check  If  pivot  point  already  normalized  to  1 

3  1 f ( a ( g * c ) . e a . 1 . ) go  to  5 
c  normalize  pivot  row 

alpha  =  a(g*c ) 

do  6  k  =  c«  n 

a(g*k)=a(g*k)/alph3 

6  1 f ( aos( a ( g» k ) ) . le» zero >  a(g*k)=0.0 

c  check  If  just  normalized  pivot  In  last  row 

5  1 f (g «ge .m ) return 

c  zero  the  elements  below  the  pivot 
gplus l=g*l 
do  7  p=gplusl»m 
b  =  a( p  •  c  ) 

1f(abs(a(p*c)).le.zero>  a(D«c>=0.0 
1 f < abs ( a (p « c ) ) .eq.O .0  )  go  to  7 
do  8  k=c»n 

8  a  (  p*  k  )  b*  a<  g*  k  )  «■  a(  p»  k ) 

7  continue 
1f(g.ge«m)  return 
g  =  q  +  l 

go  to  l 
end 
c 

c - - - - - 

c 

subroutine  pr1nt(anscol*ansrow) 
c 

c  subroutine  print  prints  the  entire  ans  matrix 
c 

c - - - - - - - - - - 

c 

c 

c  prints  ansrow  rows  by  anscol  columns 
Integer  a (40 *90 >» ans col » ansrow 
common  a *ans ( 90  *  1 90 ) 

1 1 1  =  1 

1  1t2=anscol 

4 f ( ( 1  t2-1tl >.gt.9  )  go  to  2 
1/(1t2.eq.1tl)  return 
c  less  than  10  columns  left  to  print 
wr1te(S*500) 
do  3  1=l*ansrow 

3  writ e< 6* 501 )  < ans ( 1 ♦ J ) » f = 1 t 1 » 1 t 2 > 


return 

2  1t2=1tl*9 

c  more  thrn  10  columns  left  to  print 
wr 1 te (6  *500  ) 
do  4  1=l*ansrow 

4  wr 1 te ( 6 *501 )  (ans (  1 , j >  ,  J  =  1 1 1  *  1 1 2 > 

1tl=1t2*l 
go  to  1 

500  f ormat ( 1 x ) 

501  format ( 1 x ♦ 1 0 ( e 1 1 . 4 , 1 x  )  ) 
end 

c 

- - 

subroutine  prntl(hdr,acll»acl2*arwl,arw2) 
c 

c  subroutine  prntl  prints  only  the  desired  part  of  the  ans  matrix 
c  describing  the  port  equations  along  with  the  column  headlnqs 
c 

c - - - - - 

c 

Integer  a(40*90)*header(320)«acll*acl2«arMl*arw2*hdr 
common  a *ans (90  *  180 ) *  header 
1tm2=acll-l 
1 1 1  =  1 

1  1 1 2  =  hdr 

1  f <  < 1 1 2  —  1 1 1 ) • gt • 1 9)  go  to  2 
1 f ( 1 1 2 •eq.i 1 1 )  return 
c  less  or  equal  10  columns  to  print 

wr1te(6,500>  (header(i)*1=1tl»1t2i 
1 tml=1 tm2*l 
do  3  1=arwl»arw2 

3  wr1te(6»501)  (ans ( 1 « j ) * j =1 tml « ac 12 ) 
return 

2  1 1 2=  1 1 1>  19 

c  more  than  10  columns  to  print 

wrl te( 6*500 )( header (1 ) ,1=1 tl«1t 2) 

1 tml=1tm2*l 
1tm2=1tml*9 
do  4  1=arwl»arw2 

4  wr 1 te(6  *501 )  ( ans ( 1  * j ) * j  =  1 1 ml  *  1 1 m2 ) 

1 1 1-1 1 2  + 1 

go  to  1 

500  formate lh0*10(4x*al» 12*5x>> 

501  format(lh0,10(ell.4»lx) ) 
end 


#eor 


oo  ^  O' 


c 

c 

c 

c 

c 

c 

c 

c 
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Subroutine  Lemke  for  solving  the  complementarity  problem. 
It  utilizes  the  subroutines  (matrix*  Inltla*  newbas*  sort, 
pivot  and  pprint. 


subroutine  temke(n) 


common/004/am»q*w*z»ll*b»nli*nl2*a*nel*ne2*1r*mbas1s*1zr 
dimension  am(50*50>*q(50>«b<50*50)*a(50>*mbas1s<10Q> 
dimension  w(50)*zC50) 
description  of  parameters  In  common 

am  a  two  dimensional  array  containing  the 

elements  of  matrix  m. 

q  a  singly  subscripted  array  containing  the 

elements  of  vector  q. 

11  an  Inteqer  variable  Indicating  the  number  of 

Iterations  taken  for  each  problem, 
b  a  two  dimensional  array  containing  the 

elements  of  the  Inverse  of  the  current  basis, 
w  a  singly  subscripted  array  containing  the  values 

of  w  variables  In  each  solution, 
z  a  singly  subscripted  array  containing  the  values 

of  z  variables  In  each  solution, 
nil  an  Integer  variable  taking  value  1  or  2  depend¬ 

ing  on  whether  variable  w  or  z  leaves  the  basis 
nel  similar  to  nil  but  Indicates  variable  entering 

nl2  an  Integer  variable  Indicating  what  component 

of  w  or  z  variable  leaves  the  basis. 
ne2  similar  to  nl2  but  Indicates  variable  entering 

a  a  singly  subscripted  array  containing  the 

elements  of  the  transformed  column  that  Is 
entering  the  basis. 

1r  an  Integer  variable  denoting  the  pivot  row  at 

each  Iteration,  also  used  to  Indicate  termina¬ 
tion  of  a  problem  by  giving  it  a  value  of  1000. 
mbasls  a  singly  subscripted  array-lndl cator  for  the 
basic  variables,  two  Indicators  are  used  for 
each  oaslc  variable-one  Indicating  whether 
It  Is  a  w  or  z  and  another  Indicating  what 
component  of  w  or  z. 

n  Integer  variable  Indicating  problem  size 

Initialize  basis  Inverse, 
do  9  J  =  l»n 
do  7  1  =  1  *n 
1 f ( 1 . eq. 3  >  QO  to  8 
b< 1 ♦ J  >  =0 .0 
go  t  o  7 
b(  1 »  3 )  =  1 .0 
cont Inue 
cont Inue 

parameter  n  Indicates  the  problem  size 
cal  l  In  It 1a(n) 

since  for  any  problem  termination  can  occur  In  Inltla* 
newbas  or  sort  subrout  1 ne *  the  value  of  1r  Is  matched  with 
1000  to  check  whether  to  continue  or  go  to  next  problem. 

1  f ( 1  r. eq. 1  0 00  )  return 


53 


50  call  neMbas(n) 

1f(1r.eq#1000)return 
call  sort(n) 

1f< 1r.ge.l000)return 
call  plvot(n) 
go  to  50 
end 


suoroutlne  Inltla(n) 

c  purpose-to  find  the  Initial  almost  complementar y  solution 
c  by  adding  an  artificial  variable  zo. 

c 

c - - - - - - - - 

c 

common/ 0  04/ am  to «w  «  z  * 1 1 *b  »nl 1 *nl 2*a  *ne l»ne2  tl r * mbasl s  1 1 
dimension  am(50»50)»q(50)»b(50»50)»a(50)tmbas1s(100> 
dimension  w(50)«z(50) 
c  set  zO  equal  to  the  most  negative  q(1> 

1=1 

3=2 

9  IfCqCI)  *le.  q(j))go  to  18 

1  =  3 

18  J  =  J ♦ 1 

1  f  <  j  •  le.  n)go  to  n 
c  update  q  vector 
1  r=  1 

tl=-q<1r> 

1  f ( 1 1 •  le.O . 0 )  qo  to  1000 
do  10  1  =  1  »n 
q(  1  )  =  q(  1>«-tl 
10  continue 
q  (  1  r )  =t  1 

c  update  basis  Inverse  and  Indicator  vector 
c  of  basic  variables, 
do  12  J  =  1 » n 
b( 3  * 1r) =-1*0 
w<3  >=q( J ) 
z ( J )  =  0 .0 
mbasl s( 3  >  =  1 

I  =n*  j 

mbas  Is ( l )  =  3 
12  continue 
1 zr =1r 
n  11  =  1 
l=n*1 r 
nl2=1r 

mbasl s( 1 r) =3 
mbasl s ( l >=0 
w( 1 r ) =0. 0 
z0  =  q(  1r ) 

II  =  1 
return 

c 

c 

1000  1 r=1000 

do  1010  1=1 *n 
mbasl s( 1 ) =1 
mbas 1  s  < 1*n  >  =1 
1010  w(1)=q(1) 


return 


end 


subroutine  newbas(n) 
c 

c  purpose  -  to  find  the  new  basis  column  to  enter  In 
c  terms  of  the  current  basis* 

c - — - - - - - - 

c 

common /0  04/amt  q* w « ?  t 1 1 t b  t  nil  *  nl2»  a*  nel *  ne2 1 1 r t  mbas 1 s 
Dimension  am<50t5O>tO<5O)tb<5Ot5O)ta(5O)tmbas1s(lOO) 
dimension  w(50)tz<50> 

c  If  nil  Is  neither  1  nor  2  then  the  variable  20  leaves  the 
c  basis  Indicating  termination  with  a  complementary  solution 
1f(nll  *eq*  l)go  to  20 
1f(nll  *ea*2)go  to  21 

c  If  the  complementary  solution  and  the  number  of  Iterations 
c  are  to  be  printed  set  1pp=l  In  the  following  statement 
c 

1pp  =  0 

If(lpp.eq.l)  then 
call  pprlnt(n) 
endlf 
1 r=1000 
return 

20  nel=2 
ne2=nl2 

c  update  new  basic  column  by  multiplying  by  basis  Inverse, 
do  26  1 =1 »n 
tl=0.0 
do  28  J  =  1  »n 

28  tl=tl-b< It J)*am< Jtne2> 
a  <  1  >  =t  1 

26  continue 
return 

21  ne 1=1 
ne2=nl2 

do  29  1  =  1 1  n 
a(1 )=b(1 tne2> 

29  continue 
return 
end 

c 

c---- - - - - - - - - - - 

c 

subroutine  sort(n) 

c  purpose  -  to  find  the  pivot  row  for  next  Iteration  by  the 
c  use  of  (simplex-type)  minimum  ratio  rule* 

c 

C---------e---------------------- ----------- ------------------ 

C 

Common/004/amtatwt2*llt0»nlltnl2tatneltne2t1rtmbas1st1zr 
dimension  am(50t50)tq(50)tb(50«50)ta(50>tmbas1s<100) 
dimension  w(50)tz(50) 
amax  =  abs( a  1 1 ) » 
do  101=2tn 

1 f ( amax . ge. abs ( a ( 1 > > ) go  to  10 
amax=abs(a<  1 ) ) 


10  continue 

c  set  tol =  ama k *2* 0* * ( - (b-1 1 ) )  where  b  Is  the  number  of 
c  bits  In  the  floating  point  mantissa  as  clasen  suggests* 
tol=amax*2.0**<-27) 

1  =  1 

52  1 f ( a(  1) *gt*  tol) go  to  51 
1  =  1*1 

1 f ( 1 *g t  *n ) go  to  9 
go  to  52 
51  t l  =  q  < 1  )  /a ( 1  ) 

1  r  =  1 

55  1=1*1 

1 f  < 1 *gt  *n  )go  to  56 
If (a( 1 ) . gt . tol > go  to  54 
go  to  55 
54  t2=g( 1 )/a( 1 ) 

1 f C  t  2. ge. 1 1 ) go  to  55 

I  r  =  1 

I I  =  t2 

go  to  55 

56  return 

9  1 f ( q( 1 2 r ) . gt . t ol ) go  to  57 
call  pprlnt(n) 

1 r=1000 
return 

c  failure  of  the  ratio  rule  Indicates  termination  with 
c  no  complementary  solution* 

57  print  250 

250  f ormat < 5x. 37hproblem  has  no  complementary  solution) 
print  251.11 

251  f ormat( lOx. 13h1terat Ion  no*. 14) 

1r =1001 

return 

end 

c 


c 

subroutine  plvot(n) 
c 

c  purpose  -  to  perform  the  pivot  operation  by  updating  the 
c  Inverse  of  the  basis  and  q  vector* 

c 


common/004/ am. q»w.2«ll»b. nil. nl2. a. nel»ne2.1r. mb asls 
dimension  am(50«50).q(50)*b(50.50).a(50)*irbas1s(100) 
dimension  w(50).z(50) 
do  30  1=1. n 

30  b(1r.1)=b(1r«1)/a(1r) 
q( 1r) =q( 1r )/a( 1 r) 

do  31  1 =l»n 
1 f ( 1 *eq . 1 r )qo  to  31 
q( 1)  =  q( 1 ) - q  < 1r)*a( 1) 
do  32  j  =1  .n 

b(1»J)=b(1»J)-b(1r»j)*a(1> 

32  continue 

31  continue 

c  update  the  Indicator  vector  of  basic  variables 
n 1 1 =mbas 1 s  < 1 r ) 
l  =  n*  1  r 

nl2  =  moas1s(  l) 


mbasl s ( 1 r )=nel 
mbasl s< l ) =  ne2 
11=11*1 
return 
end 


subroutine  pprlnt(n) 
c  purpose  -  to  print  the  current  solution  to  complementary 
c  problem  and  the  Iteration  number. 

c 

c- - - - - - - - - - 

c 

common/004/am»c*w*z«ll»b*nll*nl2«a»nel*ne2*1r*mbas1s 
dimension  am<50*50)«q(50)»b(50»50)*a(5Ci>*mbas1s(100> 
dimension  w(50)«z(50) 
c  zero  the  variable  values, 
do  35  1  =  1  *  n 
w ( 1 )  =0 • 0 
35  z(1)=0.0 


1  =  n*l 

J  =  1 

42  kl=mbas1s<1> 
k2=mbas1s(  j  ) 

1 f ( q< J ) • ge. 0 .0 ) go  to  45 
qC  J  >  =  0.0 

45  If (k2.eq.l)go  to  40 
z(kl)=q( j) 
go  to  41 

40  w(kl)=q(}> 

41  1=1*1 

J  =  J*1 

1 f ( J  .  le  .n  )  go  to  42 
do  44  1=l»n 

write  (6*601 )  1 «w( 1 ) » z(  1  ) 

601  f orma t ( 1 x« 1 1 « 2f 14. 4  ) 

44  continue 

write (6*602)11 

602  format( lx* 'number  of  pivots  required'*  13) 
return 


c - - 

c  subroutine  ambset  for  quantizing  the  test  node  voltaqes  into 
c  ranges  centered  around  the  voltage  value  due  to  some  fault 
c  condition* 

c 

c  Input  and  Internal  variables 

c 

c  nvt  number  of  test  nodes  -  Input 

c  nft  number  of  faults  -  Input 

c  nff  vector  containing  fault  numbers  as  sbeclfled  In  the  proqram 

c  Input 

c  ch  array  containing  branch  names  (4  characters  each) 

c  n4  starting  location  of  the  test  nodes  In  the  array  ch 

c  amr  voltage  range  -  Input*  default=lv 

c  rng  array  containing  boundaries  of  voltage  ranges 

c  lout  array  used  In  printing  the  one-zero  form  of  the 

c  ambiguity  set-fault  table* 

c  nset  no  of  sets  In  the  corresponding  test  nodes* 

c  ans  array  containing  the  vt-fault  table* 

c  aa  Integer  working  array 

c  am  real  worklnq  array 

c 

c- - - - - - - - - - - - 

subroutine  ambset(nvt*nft»ch»n4*amr*nff> 

dimension  rng(2*20)« 1out(20*40)»nset(40)*1cntr(20> 

dimension  nff<20) 

common  aa( 40 * 90 ) « ans ( 90 « 180 ) 

common/004/ am (50*50) 

Integer  aa«ch(2*90) 
c 

If (amr.eq.0.)  amr=l. 
zero=10.**<-8. ) 

c  clear  flags  and  ambiguity  set  table  taa  matrix) 
do  20  1=1*40 
do  20  5=1*40 
aa( 1  *  J ) =0 
20  continue 

c  set  fault  binary  codes  In  aa  column  40 
aa(l*40)=l 
do  5  J  =  2 » nf t 
aa< J»40)=aa( J-l»40)*2 
5  continue 

c  scan  nodes  In  v-f  table  (ans  matrix) 
do  200  n=l*nvt 

c  generate  table  of  differences  In  am  matrix 
do  10  J  =l*nf t 
k 1=  J  *1 

do  9  k=kl«nft 
am(J*k)=ans(n*J)-ans(n*k) 

9  am(k*])=am(j*k) 

10  continue 
c  clear  output  matrix 
do  25  1=1*20 
do  25  J=1 »nf t 
25  1out(1«J)=0 

c  set  ambiguity  set  1=0  for  node  n 
1  =  0 

c  scan  faults  vertically 
do  100  J  =  l,nf t 

c  if  fault  1  has  been  scanned*  skip  and  search 


e  for  a  new  center  to  the  next  ambiguity  set 
1 f < aa< J * 39 ) «ne« 0 > go  to  100 
c  update  flags  and  amolgulty  set  Index 
aa  <  j  1 39) =1 
aa< J t 38) = j 
UU1 

aa( J*37)=l 

c  set  Initial  amblg  set  code  and  outDut  matrix 
aa( l«n)=aa( 3*40) 
lout ( I* ] )-l 

c  define  amolg  set  range  and  center 
rng(l»l)=ans(ntj ) -am r/2 
rng(2»l)=ans(n»j )  ♦a mr/2 
1cntr(l)=j 

c  scan  faults  horizontally 
k  1  =  3  ♦  1 

c  reset  flags  of  updating  overlapped  ranges 
If lagl=0 
If lag2=0 
do  90  k  =  1 »nf t 

c  If  fault  k  out  of  range  skip  It 
1  f ( k.eq.  J )  go  to  90 

1 f ( abs( am( j * k) ) • ge« ( amr/2) )  go  to  90 
c  If  fault  k  has  been  scanned  check  If  It 
c  Is  closer  to  the  old  or  the  new  center 
1ct=aa(k*38) 

1 s t  =aa( k  *37) 

1 f { aa( k» 39) .ne« 1 ) go  to  60 
w=abs(am(j«k))-abs(am(1ct«k)> 

If ( w.ge.zero)  go  to  50 
c  If  fault  Is  closer  to  the  new  center  add 
c  It  to  the  new  set  ano  delete  It  from  the 
c  old  one  then  update  flags  and  parameters 
aa(l*n)=aa(l*n>.or.aa(k«40) 
a  at  1 s t  *  n)  =  aa( 1st*n)-aa(k*40) 
aa  Ck « 38 )  =  ) 
aa( k«  37)  =  l 
1  ou  t  <  l  *  k  )  - 1 

c  update  amblg  set  ranges 

50  1 f  (  am( 1 ct i ] ) • le« zero ) qo  to  55 
If  <1  flagl.ecul )  go  to  90 
c  set  flag  to  acknowledge  updating  positive 
c  side  overlap 
1 f lagl=l 

rng(2* l )=ans(n*  J )eabs(am(J*1ct))/2. 
rng(l*1st)= an s(n*1ct)-abs( am (J»1ct>)/2. 
go  to  90 

55  If (1 f lag2.eg.l >go  to  90 
If lag2=l 

rng<ltl)= an s(n*J)-abs(am(1ct*J))/2« 
rng( 2 f 1 s t  >  =ans ( n ♦ 1 ct ) ♦abs ( am< 1 c t » J  > ) /2. 
qo  to  90 
c 

e  If  fault  has  not  been  scanned  and  within  range 
c  add  It  to  the  current  set  and  update  parameters 
60  aa  1 1 »n  )  =  aa ( l »n )  .or «aa ( k % 40 ) 

1 Out< l» k)=l 
aa ( k  •  39)  =  1 
aalkf 38)=J 
aa( k  t  37) =1 


90  continue 

1  f  ( l  »eq  •  1  >  go  to  100 
do  95  11=1.1 

1f((rng(l»l).tt.rng(2»tl)).and.(rng(l»L).gt.rng(lfll)))  then 
rnq<  1 . 1)  =(  rnq(  1 «  l)  ♦rng(  2.  ID  >  /  ?  • 
rnq ( 2 » 1 1 ) =rnq (  1 « l ) 

else1f<(rng(2.l).gt«rng(l»ll)>.and*(rng(2.l)*lt«rng(2.ll))>  then 

rng(2»l)=(rng(2*l>*rng<l»tl>)/2. 

rng< 1 . I l ) =rng( 2 . 1 ) 

endl  f 

95  continue 
100  continue 
nset  (  rt)  =  l 

c  write  out  centers,  ranges  and  amlqulty  set  codes 
c  for  node  n 

wr1te<6»220>  ch{l.n*n4).ch(2»n*n4) 
wr1te(6.230) 
do  110  1  =  1.  I 
«  =  1cntr(1 ) 

1 cnt  r< 1 ) =nf  f <  m) 

110  wr1te«6.240)  1,1  cntr(1).<rngCj,1),j  =  i»2), (lout (1 .J >» j  =  l.nf  t  > 

c  clear  flags 

do  115  1=1. nft 
do  115  1=37.39 
115  aa( 1 . j ) =0 
200  continue 

220  f ormat < //2x » 'node  ->  1  ’»2a2) 

230  format(//2x**set'« '*»»center*»3x.'from*»10x**to'.15x»'set  code') 
240  format</14,l6.2(2x.el0.3>»5x.40(12>> 

call  ftcode (nvt. nft. nset* ch.n4»nff) 
return 
end 


#eor 


c  subroutine  ftcode  tor  generating  a  numeric  code  fault  dlctlona- 

c  ry  based  on  the  voltage  ranges  provided  by  the  subroutine 

c  "ambset". 
c 

c  Input  and  Internal  variables 

c 

c  nvt  number  of  test  nodes  -  Input 

c  nft  number  of  faults  -  Input 

c  nset  number  of  ranges  In  the  corresponding  test  node  -  incut 

c  ch  array  containing  branch  names  -  Input 

c  n4  starting  location  of  test  nodes  In  ch  -  Input 

c  nff  vector  containing  fault  numbers  -  Input 

c  Icod  array  containing  the  numeric  codes  on  output 

c  aa  Integer  working  array 

c 


subroutine  ftcode(nvt»nfttnset*chtn4»nff) 
dimension  nset(40)«1ax(40)«1cod(40*40)t1temp(40) 
Integer  ch(2»90) 

Integer  a«b«nff(30) 
common  a<40*90) 
c 
c 

c  set  node  Index  to  zero  and  find  the  node  with  maximum 
c  number  of  ambiguity  sets 
1 debug=0 
n=0 

lmax=nset (1  ) 
do  10  1=l»nvt 

10  lmax  =  max0( lmaxfnset( 1 ) ) 
lmaxl=lmax 

1f< Idebug.eq.O)  go  to  201 
wr1te(6«ll)  Imaxl 

11  format(2x»» lmaxl=*t120) 

c  move  the  cor  responding  set  codes  In  ’a*  to  temporary 
e  storage  and  Identify  the  node 

201  do  20  1= 1 »nvt 

1 f (nset( 1 ) . eq. Imax)  go  to  1R 

20  continue 
18  11=1 

nsetl 11 >  =  0 
do  19  J=l«lmaxl 
19  ItempC J)=a()*11) 
n=n+ 1 

1 f ( Idebug.ea. 0  )  go  to  202 
wr1te(8*41)  n 
wr1te(6» 21)  11 

21  format(2x» »node» » 15) 

202  nn=n 

c  Initialize  the  fault  code  matrix  1cod(*«») 
c  by  the  ambiguity  set  Indeces  of  the  first  node 
do  22  1=l»lmaxl 

22  1cod(1tl)=1 
lf=lmaxl 

If(nvt.gt.l)  go  to  30 
do  23  j=l»lmaxl 

23  1  ax ( j )  =  1  temp ( j  ) 

=  *  ma  x  1 


go  to  ICO 

c  find  the  next  maximum 
30  Uai:n$et(l) 
do  35  1=l»nvt 

35  l  ma  x  =ma  x  0  <  l  ma  x  t  n  se  t  {  1  )  ) 
lmax2= Imax 

1 f ( ldebug.eq.0  )  go  to  203 
wr1te<6»36>  lmax2 

36  fo.-mat(2x.»lmax2=*t120> 

203  do  40  1=l»nvt 

1 f Inset <<). eq« Imax)  go  to  42 

40  continue 

42  nsec<1)=0 

12  =  1 
n  =  n*l 
nn=nn»l 

1 f ( 1 deoug .eg. 0 )  go  to  204 
mH  te  (6*41 )  n 
wr1te<6f21>  12 

41  format(2x»,n=,#l4) 

c  find  the  Intersection  of  the  am  sets  by  logically 
c  andlng  the  set  codes 

204  1=1 

do  50  j=l?lmaxl 

c  set  a  flag  to  detect  first  Intersection 
If  1  =  0 

do  50  k=l»lmax2 
1jx(l)=1temp(j).and.a(kt12) 

IMIax(l).eq.O)  go  to  50 
c  otherwise  update  l  and  If  and  Insert  the  code 
c  1 cod< J » m) « m= 1 • n-1  In  the  next  locations  (J»l) 
c  pushing  own  the  rest  of  the  codes  up  to  the 
c  location  if.  then  add  the  set  Index  k  In  the 
c  location  1cod(]*n)  corresponding  to  set  J  -node  n 
J l=lf-l*2 
nl =n-l 

c  skip  pushing  codes  after  first  Intersection  only 
IfCIfl.ne.O)  go  to  44 
1cod(l»n)=k 
If  1  =  1 
go  to  47 

44  do  45  11=1*J1 
do  45  m=l»nl 

45  1cod<lf-l1*2»«)=1cod(lf-l1*lfn) 

Icodt It n) =k 

If =lf *1 

47  1=1*1 

1f( l. le.nft)  go  to  50 
1  =  1-1 
go  to  100 
50  continue 
1  =  1-1 

1 f ( Idebug. eq.0)  go  to  205 
do  48  ll=lfl 

48  wr1te(6f49)  < 1 c rd ( l l * kk > »kk =1 *n > 

49  format(2xf *code’ il015) 

c  check  If  the  new  node  Is  redundant 

205  1 f ( l .ne .ma x 0 ( Ima xl « Ima x2  )  )  go  to  70 

c  If  the  number  of  sets  rsultlng  from  Inersctton 
c  Is  not  Increased  Ignore  the  node  and  go  to  pick 
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c  a  new  node 
n=n-l 

1fCnn.eq.nvt)  qo  to  70 
go  to  30 

c  check  If  ill  faults  have  been  Isolated  or  all  nodes 
c  are  exhausted 

70  1 f ( l. eq. nf t • or. nn. eq. nv t )  go  to  100 

c  move  result  of  Intersection  to  temporary  storage 
lmaxl=l 

do  80  1=l,lmaxl 

80  1  temp  C 1 )  =  1  a  x  <  1  > 
go  to  30 

100  wr1te<6*120) 

wr1te<6*130)  Cchei,n4*m),che2»n4*m>*«=l»n) 
do  95  1=1,1 
do  82  j=l»nf t 
1tl=1ax(1).and.a(J ,40) 

If < 1tl.eo.1axC  1) )  go  to  81 
If C ltl.eq.0  )  go  to  82 
go  to  83 

81  ur1teC6, 131 >  nffCj) 
wr1tel6»140>  C 1 ccdC 1 »m ) , m=l ,n) 
go  to  95 

82  continue 

83  m=0 

do  85  J=l,nft 

1fCfloatC1axC1))/2.0.eq.floatC1axC1)/2))  go  to  84 
m=m  +  l 

a(k,41)=nffC J) 

84  1axC1)=1axCl)/2 

85  continue 

wr1teC6,150)  C a C k  ,41 )  ,k  =  1  ,m  ) 
ur1teC6,140  )  C 1  code  1 , j ) , j =1 ,n> 

95  continue 

120  f ormat C //2x »* fault ' ,23x» « fault  code*) 

130  formatC/30x,20C2a2,lx>  ) 

131  f ormat  C2x,*f»*1?) 

140  formatC30x,20C  13,2x) ) 

150  formatC2x,»f « , 10 C 1 2, » , » ) ) 
stop 
end 


»eor 


Appendix  S 
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MIFTCOD:  A  Subroutine  for  Fault  Isolation  Under  Multiple  Test  Input  Conditions. 

In  case  of  applying  different  test  input  conditions  to  improve  the  level  of  isolation, 
there  are  two  ways  of  handling  the  fault  dictionary.  The  first  way  is  to  produce  a  fault 
code  combining  the  input  conditions,  and  that  is  to  be  used  if  the  measurement  data  is 
available  all  at  once.  The  other  way  is  to  generate  a  dictionary  for  every  input.  Then 
the  isolation  is  achieved  by  correlating  the  dictionaries  of  the  different  inputs. 

The  subroutine  MIFTCOD  accepts  the  result  o»  simulation  for  all  input  conditions 
after  being  classified  into  ambiguity  ranges  and  produces  the  required  fault  codes  for 
the  individual  inputs  and  the  codes  combining  all  inputs.  This  subroutine  is  not  util¬ 
ized  by  the  program  HAFDIC.  Some  extra  effort  has  to  be  done  in  writing  a  small  pro¬ 
gram  that  takes  the  results  of  HAFDIC  simulation  under  different  input  conditions  and 
then  calls  MIFICD  to  perform  the  required  isolation.  The  program  BASTEX  shown 
next  is  an  example  of  this. 

Example  for  Using  MIFTCD: 

Table  A3-1  contains  the  ambiguity  sets  of  faults  simulated  in  a  circuit  taken  from 
(4],  where  two  test  inputs  are  used  separately.  To  be  able  to  supply  the  MIFTCD  sub¬ 
routine  with  the  ambiguity  sets  in  the  required  form,  the  ith  fault  is  assigned  the 
integer  weight  2‘  (see  BASTEX  program).  This  means  that  the  integer  word  represent¬ 
ing  this  fault  will  have  the  ith  bit  equal  to  one,  while  all  other  bits  are  zero.  Adding  or 
logical  ORing  of  two  faults  will  simply  mean  that  there  will  be  two  bits  equal  to  one  in 
the  integer  word,  and  so  on  for  more  faults.  The  program  output  is  shown  next  to  the 
ambiguity  sets  table  followed  by  the  listing  of  MIFTCD. 


4 


4 


rUSTEX 


c  This  program  arranges  the  ambiguity  sets  of  faults  shown  in 

c  the  accompanying  table  in  a  two  dimensional  array  that  orovld 

c  the  required  input  to  the  subroutine  MIFTCO. 

c 

program  bas t e x ( input , out  out , t ape5  =  i nput ♦ tape6=  out pu t ) 
c  this  program  generates  the  ambig  sets  of  bastian7s  oaoer 
c  example  in  the  matrix  a* 

dimension  nse t < 4 » 20 ) , 1 vt ( 20 ) 

Integer  a,c<20) 
common  a(50»100) 
c( 1) =1 

do  10  1=2*20 
10  c( 1 ) =2* c( 1-1 ) 

1 cs=c ( 1 > 
do  5  1=2*20 
5  1cs  =  1cs  +  c ( 1 ) 

a(l*l)=c(3)*e (6)*c(7)*c(19) 

a  <  2  » 1 )  =  1 cs-a( 1*1) 

a  < 1 ,2  >  =c<  3) 

a (2*2)=c (7) 

a( 3, 2  >  =  c ( 15 > 

a(4*2)=c<16> 

a<5*2)=1cs-a<l»2)-a<2»2>-a(3»2>-a<4*2) 

a(l*3)=c(3)+c(6)-»c(7)*>c(15) 

a  ( 2  » 3 )  =  1 cs-a ( 1 , 3 ) 

a(l*4)=a(l,3) 

a<2*4) =c< 19) 

a (3*4)=c (20 ) 

a(4*4)=1cs-a(l*4)-a(2*4)-a(3,4) 

a(l*5)=a(l*3)+c(19) 

a ( 2*5) = 1 cs-a (1*5) 

a( 1*6) =a( 1*5) 

a(2»6)=a(2*5) 

a(l, 7)  =a (1,5) 

a(2*7)=a(2*5) 

a(l*8)=a(l,3)+e(10)+c(12)+c(l9) 

a( 2  *8 )  =  1  cs-a( 1*  8) 

a( 1,9) =c(2) 

a<2 , 9 ) =c (5) 

a  <  3  *9  )  =  c  <  1 3  ) 

a(4*9)=c(14) 

a(5,9)=1cs-a(l,9)-a(2,9)-a(3,9)-a(4,9) 

a(l*10)=c(2)+c(4)+c(5)+c(13) 

a(2,l0)=1cs-a( 1*10) 

a(l,ll)=a(l»10) 

a(2,ll)=c(9) 

a(3,ll)=c(17) 

a( 4,  ID  sc  (13) 

a(5,il)=1cs-a(l*ll)-a(2»ll)-a(3»ll)-a(4,ll) 

a(l«12)=a<l*10)*c(8)*-c<9)*-c(17> 

a(2«12)=1cs-a(i,l2) 

a( 1 , 13)  =  a( 1  *  12) 

a( 2* 13) =a( 2*12) 

a(l,l4)=a(l«12) 

a( 2* 14) =a (2  *  12 ) 

a(  1 , 15) =a( 1 , 12 ) 

a (2*15 )  =  a(2  *12 ) 

a( 1  *  16 ) =a( 1  *  12 ) 


a  ( 2 *  16 )  =e  <  1 1 ) 

a(3»16)=1cs-a(l»16)-a(2«16> 

1 vt ( 1 ) =1 l 
1vt<2>=8 
1 vt ( 3) =5 
Ivt (4)=2 
Ivt <5>=27 
Ivt ( 6 ) =26 
1vt<7>=33 
Ivt ( 8) =16 
nsetc 1 ♦ 1 ) =2 
nset (1» 2>=5 
nset( It 3)=2 
nset  < If  4) =4 
nset< 1 f  5  >  =  2 
nset ( If  6) =2 
nset ( lf7)=2 
nsetdf  8)=2 
nset  < 2t 1 >=5 
nset (2t 2)=2 
nset<2»3)=5 
nset(2«4>=2 
nse t ( 2f 5>  =2 
nset<  2f6)  =  2 
nset(2f7)=2 
nset(2f8)=3 
nvt  =8 
nf t=20 
n  1  =  2 

call  mlftcd(nvtflvtfnlfnsetfnftfC) 

stop 

end 


3,6,7,15,19 

2 


3 

2,4,5,13 


3,8,7, 

2,4,5, 


3,6,7,15 

2,4,5,8,9,13,17 


3,6,7,15,19 

2,4,5,8,9,13,17 


3,6,7,15,19 

2,4,5,8,9,13,17 


3,6,7,15,19 

2,4,5,8,9,15,17 


nominal 

nominal 


nominal 

nominal 


nominal 

5 


7 

nominal 


nominal 


nominal 


nominal 

9 

17 

18 

nominal 

19 

nominal 


nominal 

nominal 


nominal 

nominal 


nominal 

nominal 


nominal 


3,6,7,10,12,15,19  nominal 
2,4,5,8,9,15,17  11 


nominal 


combined  Inputs 


f  3, 
f  7  * 
f  6* 
f  19, 
f  2, 
f  5  * 
f  13, 
f  14, 
f  15, 
f  16, 
f  4, 
f  9, 
fl7, 
fia, 
f  20, 
f  8, 
f 10,12, 
till 
f  1, 

Inputl 
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vll«lvll,2v  8,lv  8»2v  5»lv  5»2v  2,lv  2 »  2v  1 6 » 1  v  1 6 , 2  v 


15  12  1 
1  5  2  2  1 
1  5  5  2  1 

1  5  5  2  2 

2  15  12 
2  2  5  12 
2  3  5  1  2 
2  4  5  2  2 
2  5  3  2  1 
2  5  4  2  2 
2  5  5  1  2 
2  5  5  2  2 
2  5  5  2  2 
2  5  5  2  2 
2  5  5  2  2 
2  5  5  2  2 
2  5  5  2  2 
2  5  5  2  2 
2  5  5  2  2 

separate  Input  cedes 


5  12  1" 

5  12  15 

5  12  15 

5  2  2  1  ? 

14  12  1 

14  12! 

14  12  1 

5  4  2  2  5 

5  1  2  1 

5  4  2  2  7 

14  1  2  1 

2  4  12  1 

3  4  1.2  1 

4  4  2  2  5 

5  3  2  2  7 

5  4  12  1 

5  4  2  1  3 

5  4  2  2  2 

5  4  2  2  3 


fault  code 

fault 
f  3, 

f  7, 

f  6, 

f  1 9, 


vll  v  8  v  5  v  2  vl6  v 
11111 
12  111 

15  111 

1  5  2  2  1 


subroutine  mlftcd  for  generating  a  numeric  code  fault  d 
under  multiple  test  Input  conditions* 

Input  data 


the  ambiguity  sets  of  faults  have  to  be  stored  In  a  two 
dimensional  Integer  array  such  that  every  fault  Is  represent ed 
by  a  single  bit  In  an  Integer  word.  Every  ambiguity  set  of 
faults  Is  then  contained  In  one  Integer  word.  Therefore  the 
number  of  faults  Is  limited  to  the  computer  word  lenqth. 

If  the  program  utilizing  mlftcd  Is  to  be  written  In  standard 
fortran  where  bit  manipulation  Is  not  available,  the  Individual 
bits  can  still  be  affected  (see  the  example  In  appendix 
fault  diagnosis  of  nonlinear  analog  circuits  vol5). 


the  arrange 

ment 

at  the 

ambiguity 

sets  of  nvt 

test 

nod^s  1  n 

the  array  a 

should  be 

as 

fol lows : 

8(1.1) 

1th 

ambl  a 

set 

o  f 

the 

first  test 

node 

-  Input 

1 

a(1 .2) 

1th 

amb  1  q 

set 

of 

t  he 

second  test 

node 

-  Input 

1 

a ( 1  .nvt ) 

1th 

amblg 

set 

of 

t  he 

last  test  node  - 

Input  1 

a  < 1 »nvt  *1 ) 

1th 

amblg 

set 

of 

the 

first  test 

node 

-  Input 

at  1 »nvt*2 ) 

1th 

amblg 

set 

o  f 

t  he 

first  test 

node 

-  Input 

2 

a( 1 . 2*nvt ) 

1th 

ambl  g 

set 

of 

the 

last  test  node  - 

input  2 

at  1 . 2*nvt ♦! ) 1 1  h 

amblg 

set 

of 

the 

first  test 

node 

-  input 

3 

. .etc  • 

input  variables 

nvt  number  of  test  nodes 

Ivt  array  containing  the  actual  numbers  of  the  test  nodes 

(for  printing  the  fault  dictionary) 
nl  number  of  Inputs 

nset  array  containing  the  numbers  of  amblg  sets  In  the 
corresponding  test  nodes 
nft  number  of  faults 

c  Integer  array  containing  Integers  ( 1 . 2 .4 « • . ?**nf t ) 


subroutine  m1ftcdtnvt.1vt.nl. nset.nft.c) 
dimension  It emD<40)»1ax<40)»1vt<nvt)*1cod (40.40) 
Integer  a. c (nf t ) . cl d( 40 . 40) 
dimension  nset(4.20).mset(20).node(20) 
common  a(50.100) 
c 

c  find  the  Intersections  of  the  ambiguity  sets  of  every 
c  Individual  node  Ignorlnq  null  Intersections 
c 

nvt  2=ni *nvt 
c  clear  lax 

do  10  1=1.40 
10  1ax( 1 >=0 

do  70  j=l»nvt 
J2=n1*<  J-l) 

1  =  1 

nl=nset< 1. J ) 
n2=nset (2. ) ) 
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do  20  1 =1 *nl 
do  20  k=l*n2 

1ax(l)=a(1»J)»anc.a<k»J*nvt) 

1 f  ( 1  ax ( l > *eq . 0  )  go  to  20 
c1d<  l*  J2«-l>  =  1 
c  1  d  ( l  •  j  2+2 )  =  k 
1=1*1 

20  continue 
1  =  1-1 

c  store  the  result  ( i ax ( i ) *  1 =1  * l )  in  col  #(nvt2*j>  for 
c  further  intersection  with  the  sets  of  the  next  Input 
do  30  1=1*1 

30  a  ( 1  »n vt 2* j ) =1  a x ( 1 ) 
inset  <  J  >  =  l 

c  Intersect  a( 1»nvt2* J)  with  next  input  a(k« j*nvt 1) 
c  however  if  no  of  Inputs  «lt*  3  skip  and  go  to  next  node 
1f(n1.lt.3>  go  to  70 
do  60  n=3*ni 
nl=nset (n*  J ) 
nvt l=nvt  *  <n-l ) 
lf  =  l 
fl!=l 

do  40  1=1*1 
1  f  l  =  0 

do  40  k=l*nl 

1ax(m)=a(kfj>nvtl)*and*a(1«j*nvt2) 

1 f ( la x (m  )  .eq.O  )  go  to  40 

Jl=lf-l*2 

nn 1 -n- 1 

if(ifl.ne.O)  go  to  34 
e  1  d<  m»n«-  J2>  =k 
1  f  1  =  1 
go  to  37 

34  do  35  11=1*31 
do  35  m1=l«nnl 

35  cid<lf-U*2*mi*j2)  =  cidC  lf-l1*l*mWJ2> 
c 1 d( m  * n*  J2)  =  k 

lf=lf*l 
37  m  =  m  +  l 

1 f ( m.gt *nf t ) go  to  90 
40  continue 

c  adjust  the  resulting  no  of  ambiq  sets 
l=m-l 

c  restore  the  resulting  sets  of  the  new  Intersection 
do  50  1=1*1 

50  a< 1 » J ♦nvt2) =1ax( 1 > 

60  continue 

c  Identify  the  no  of  comoined  Inputs  amb  sets  for  node  j 
•set ( J ) =1 
70  continue 
go  to  100 

c  if  current  node  Isolates  all  faults  acknowledge  it 
c  and  go  to  identify  faults 
80  wri te(6*90)  ivt(J) 

90  f ormat ( /2x * • node  * « 1 2  * *  is  enough’) 
c 

c  output  the  combined  inout  codes  for  the  single  sufficient  node 
c 

wr 1 te C 6 » 95 )  < 1 vt ( J ) » n, n  =  l ,ni ) 

95  formati 1 x* • vt* *  12* *  *  *  *  11) 


do  96  1=l»nft 
do  96  k  =  ltnf t 
1tl=1ax(1).and.c(k) 

1 f ( 1 1 l*ne • 1  a x ( 1  )  )  go  to  96 
wr1te<6»340)  < c 1 d< 1 ♦ m) » m=l «n1 )  ♦  k 
96  continue 
n=l 

node(l)=J 
go  to  200 

c  now  the  result  of  all  Inputs  for  every  node  5  Is  available 
c  In  a( 1 tnvt2+ j)  and  the  cor responding  no  of  ambiguity  sets 
c  Is  In  fflset(j).  output  mset(J).  perform  Intersections  of 
c  sets  belonglno  to  different  nodes  In  a  descending  order  of 
c  the  no  of  amplg  sets 
100  wr 1 te( 6* 105  >  ( mset< j > » J  =  1 »nvt ) 

105  format (/2x • ‘combined  Input  set s * /2x * 201 3) 

Ima x=mset ( 1 ) 
do  110  1=l»nvt 

110  Imax =max 0( Imax «mset < 1 ) ) 
lmaxl=lmax 
l f =lmaxl 

c  move  the  e  ••^ponding  sets  to  temporary  storage  and 
c  the  node 

do  130  1-l«nvt 

1 f ( mset (  1 ) . eq« l maxi )  go  to  135 
130  continue 

135  node< 11=1 
mset  <  1  >  =  0 
n  =  l 

nn=l 

11-1 

do  136  1=1* l maxi 

136  1temp<1)=a<1»11*nvt23 

c  load  Icodl.f.)  with  the  combined  Input  Indeces  of  the  first  node 
J2=n1*< 11-1) 
do  13d  1=l«lmaxl 
do  138  m=l*n1 

138  1cod<1*m)-c1d<1»ir*j2) 

c  find  the  next  max  and  do  the  same 

139  Imax-mset(l) 
do  140  1=l*nvt 

140  lmax=max0( lmax«mset< 1) ) 
lmax2=lmax 

1 f ( nn. eo* nvt )  go  to  190 
If Clmax2.eq«0>  stop 
do  150  1=ltnvt 

If (mset ( 1 )«eq«lmax2)  go  to  145 
150  continue 
145  mset  <  1 )  =  0 
nn-nn* 1 
n=n*  1 
12=1 

node(n) =1 
1  =  1 

J  2=n1 * (12-1) 

Jn=n1 *<n-l ) 
do  163  1=l»lmaxl 
If  1  =  0 

do  160  k=ltlmax2 

1axtl)=1temp(1)»and«a(kt12*nvt2) 


1 f ( i ax ( l ) «eq • 0 )  go  to  160 
J 1= l f-l *2 
n 1  =  n- 1 

If(lfl.ne.O)  go  to  154 
do  153  m=l«n1 

153  1  cod  < l  jn  )  =c  1  d  < k  »m>  J  2 ) 

1fl=l 

go  to  157 

154  do  155  U  =  i*Jl 
do  155  m1=l»Jn 

155  1  cod(  If- IW»  ml )  =  1cod(  I f-l 1+1 » ml) 
do  156  m=l»n1 

156  IcodC  I»m+Jn)=c1d(k«m+  J  2 ) 

If  =  Lf  <-1 

157  1fd.eq.nft)  go  to  190 
1  =  1  +  1 

160  continue 
1  =  1-1 

wr1te(6tl61)  l 

161  format (2x* ' l  =  *  *  15) 

1 f ( l .ne.maxOC Ima x 1  * Ima x2) )  go  to  170 
c  If  current  node  Is  recundant  Ignore  It  and  consider  a  new  node 
n=n-l 
go  to  139 

c  check  If  all  nodes  have  been  scanned 
170  1fCnn.eq.nvt)  go  to  190 

c  If  yes  go  to  find  the  fault  codes  otherwise  move  the  sets 
c  to  temporary  storage  and  continue 
lmaxl= l 

0 o  180  1=l«lmaxl 

180  ItempC 1)=1ax( 1) 
go  to  139 

c  write  combined  Input  fault  table 

190  wr1te(6.191) 

191  format(2x»'comb1ned  Inputs') 
wr1te(6«189) 

189  formate  1  x» '------ - - '  ) 

do  192  1  =  1  »n 
m=node  < 1 ) 

192  mset ( 1 ) =1vt  Cm ) 

wr1te<6»193)  CCmsetC 1)tnm#nm=l*n1)t1=l»n) 

193  format<4x»30('v'*12»'»'t11)> 
nn=n*n1 

c  nn  Is  the  total  no  o entries  In  a  row 
do  195  1=l»l 
do  182  J  =  1 1 nf t 
1 1 1  =  1  ax  1 1 ) . and.c C  J ) 

1 f(  1  tl.ee. 1ax< 1 ) )  qo  to  181 
1 f ( ltl.eq.0)  go  to  182 
go  to  183 

181  wr1te<6f34Q)  ( IcodC 1  * m) ♦ m  =  l «nn> ♦ j 
go  to  195 

182  continue 

c  find  unisolated  faults 

183  m=0 

do  185  J  =  1 »nf t 

1f(float(1ax(1))/2.0.eq.float(1ax(1)/2))  go  to  184 
m=m  + 1 

a( mi 100 )  =  J 

184  lax ( 1 ) =1ax ( 1 )/2 


185  continue 

wr1te<6»340)  (1coa(1«j)fj=ltnn)«(a<k«100)tk=l»m) 

195  continue 

c  proceed  to  find  the  fault  codes  based  on  the  selected  set 
c  of  nodes  In  node(j).  the  number  of  nodes  =  n.  the  no  of 
c  amblg  sets  Is  already  In  nset(nltj). 
c  the  procedure  Is  repeated  for  alt  Inputs. 

200  wr1te<6»201) 

201  f ormat ( /20x • « separa t e  Input  codes*) 
do  300  n11=l»n1 

nvt l=nvt  *  <n 1 1 -1 ) 
m=node ( 1 ) 
nl=nset ( nl 1 «m) 
lf=nl 

wr1te(6*213)  nvtl 

213  formatl 2 x* • nvtl  =  •»  15 ) 
wr1te(6»214)  1vt(m)«nl 

214  f ormat<2x* • vt * « 13» •  no  of  sets**13) 

c  load  Icod  with  the  Indeces  of  the  first  node 
c  and  move  amb  sets  to  temporary  storage 
do  215  1=l«nl 

215  1cod(1tl)=1 
do  220  1=l»nl 

220  1  temp ( 1 ) =a( 1 »m*nvt 1 ) 

If(n.ne.l)  go  to  225 

c  If  one  node  Is  enough  go  to  Identify  faults  directly, 
do  222  1 =l»nl 
222  lax ( 1 )  =  1  temp ( 1  ) 
go  to  270 

c  start  Inerseetlon 
225  do  260  nn=2*n 
1=1 

m=node<nn) 
n2=nset(n11«m) 
wrltel 6.214)  1vt(m)«n2 
do  250  j=l»nl 

c  set  flag  to  detect  first  Intersection 
1  f  1  =  0 

do  250  k=l t  n2 

1ax(l)=1temp(J ) . and. a < k ♦ m*nvt 1 ) 

1 f ( 1 ax( l ) • eq. 0 >  go  to  250 
c  update  Indeces  for  stack  pushing 
J 1=1 f-l *2 
nnl=nn-l 

c  skip  pushing  after  the  first  Intersection  only 
If(lfl.ne.O)  go  to  244 
1cod< l«nn)=k 
1fl  =  l 
go  to  247 

244  do  245  11=1«J1 
do  245  m1=i»nnl 

245  1eod<  lf-l1*2»m1)=1cod(lf-l1*l»m1> 

1cod(l»nn)=k 

If =lf *1 
247  1=1*1 

250  continue 

c  adjust  resulting  no  of  amb  sets 
1=1-1 

wrl  ,.e(6'161 )  l 
do  251  1=l«l 


251  wr 1 te( 6,  340  )  ( 1  codM  ,  J )  *  J  =  1  ,nn ) 
c  wove  result  of  Intersection  to  the  temporary  storage  and 
c  continue  to  consider  next  node 
do  255  1=1,1 
255  1  temp ( 1 ) =1a x ( 1 ) 

nl  =  l 
lf  =  nl 

260  continue 

c  now  1ax(.)  contains  all  sets  representing  Isolated  faults 
c  and  Icod  contains  corresponding  fault  codes*  proceed  to 
c  Identify  the  Isolated  faults  and  print  out  the  Icod  matrix 
write(6*370)  nil 
w  rl te( 6 , 380  ) 

c  write  the  actual  node  numbers 
do  275  1=1, n 
m=node(  1 ) 

275  mset  ( 1 >  =  1 vt <  m) 

270  wr1te(6,320> 

wr1te(6,330>  ( mset< 1 ) , 1=1 »n> 

c  match  sets  against  binary  representation  of  the  faults 
do  295  1=1,1 
do  282  J  =  1 «nf t 
1 1 1  =  1  a  x  < 1 ) • and*  c ( j ) 

1 f ( 1 1 l*eq* lax ( 1 ) )  qo  to  281 
1f(1tl*ec*0)  go  to  282 

c  If  they  do  not  match  or  andlng  not  equal  to  zero  then  there 
c  must  be  a  set  of  unisolated  faults 
go  to  283 

281  wr 1 te ( 6 , 340  >  ( 1  cod ( 1 • m) ,m=l , n) ,  J 
go  to  295 

282  continue 

c  find  unisolated  faults 

283  m=0 

do  285  j  =l«nf t 

1  ft f loat ( 1ax( 1 )) /2.0 *eq* float ( 1  ax ( 1 >/2> )  go  to  284 
m  =m*  1 

a(m,100)=j 

c  shift  olnary  code  once  to  the  right 

284  1ax(1)=1ax(1)/2 

285  continue 

wr 1 te (6, 340 )  (Icod (1,J),J=l,n),(a(k, 100),  k=l,m> 

295  continue 

c  continue  to  consider  next  Input 
300  continue 

320  f ormat ( //2x , • f au 1 1  code*) 

330  f ormat </2x, 20 ( 'vt', 12, lx), 'fault • ) 

34G  for mat(4x, 20(13, 2x) ) 

370  f ormat ( //2x ,* Input *, 1 1 ) 

380  format  (/2x,  — » ) 

return 
end 


